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ABSTRACT 

We investigate the abundance of galactic molecular hydrogen (H 2 ) in the “Evolution and 
Assembly of GaLaxies and their Environments” (EAGLE) cosmological hydrodynamic sim¬ 
ulations. We assign H 2 masses to gas particles in the simulations in post-processing using two 
different prescriptions that depend on the local dust-to-gas ratio and the interstellar radiation 
held. Both result in H 2 galaxy mass functions that agree well with observations in the local 
and high redshift Universe. The simulations reproduce the observed scaling relations between 
the mass of H 2 and the stellar mass, star formation rate and stellar surface density. Towards 
high redshifts, galaxies in the simulations display larger H 2 mass fractions and lower H 2 de¬ 
pletion timescales, also in good agreement with observations. The comoving mass density of 
H 2 in units of the critical density, Uhs, peaks at z ss 1.2 — 1.5, later than the predicted peak of 
the cosmic star formation rate activity, a.tziv2. This difference stems from the decrease in gas 
metallicity and increase in interstellar radiation held with redshift, both of which hamper H 2 
formation. We hnd that the cosmic H 2 budget is dominated by galaxies with > 10® Mq, 
star formation rates > 10 Mq yr“^ and stellar masses Msteiiar > 10^® Mq, which are readily 
observable in the optical and near-lR. The match between the H 2 properties of galaxies that 
emerge in the simulations and observations is remarkable, particularly since H 2 observations 
were not used to adjust parameters in EAGLE. 
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1 INTRODUCTION 


The past few years have brought impressive developments in sur¬ 
veys of molecular gas in reso lved and unresolved galaxies, lo 


cally and at high redshift (e.g. iLerov et alj|2008l: ISaintonge et all 
201 ll: iGenzelet alJbOld: lDaddiet alJl2010l: lYoung et alj I 2 OIII: 


Davis et al. 2011 : Genzel et al. 20131 : Boselli et all 2014) . Molec- 

ular hydrogen (H 2 ) is an important component of the interstel¬ 
lar medium (ISM) and, in galaxies l ike our own, it am ounts 
to a few percent of the stellar mass jPutman etalj|2012h . The 
exquisite resolution and sensitivity of studies of local galaxies 
have enabled important conclusions to be reached. For example, 
on kiloparsec scales, observations show a nearly linear correla¬ 
tion between the star form ation rate (SFR) surface density and 
the surface density of H 2 dBigiel et alj 1200 Sl : iBigiel et"^ [2OI0I : 


* E-mail: claudia.lagos@icrar.org 


iLerov et^l2008l : ISchruba et al.ll201 ill . At a global level (i.e. in¬ 
tegrated o ver the galaxy), the H 2 mass correlates ver y well with the 
total SFR dSaintonge et al.|[201 itlBoselli et al.ll2014ll . These obser¬ 
vations place new constraints on galaxy formation models. 

Molecular hydrogen is very difficult to observe in the ISM 
of galaxies because it lacks a dipole moment, making its emis¬ 
sion extremely weak at the typical temperature of the cold ISM. A 
widely-used tracer of H 2 is the carbon monoxide (hereafter ‘CO’) 
molecule, which is the se cond most abundant mo lecule after H 2 , 
and is easily excited (see ICarilli & Waited l2013l for a recent re¬ 
view on how CO has been used to trace H 2 at 2 < 6). Direct 
CO detections are available for relatively large samples in the lo¬ 
cal Universe (z < 0.1), from which it has b een possible to de - 
rive the 2 «OCO(l — 0) luminosity function dKeres et alj|2003h . 
where (1 — 0) is the lowest energy rotational transition. From this 
luminosity function, a nd adopting a Milky -Way like CO(l — 0)- 
H 2 conversion factor dBolatto et al.l 1201^ . it has been possible 
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2 Claudia del P. Lagos et al. 


to derive the H 2 mass function and the cosmic de nsity of H 2 , 
Ph-> = 4.3 ± 1.1 X 10^ fe Mf7) Mpc ~^, at z ~ 0.05 jKeres et"^ 
I2003l : l0hreschk:ow & Rawling jl2009h . 

Campaigns to obtain constraints on ph 2 at redshifts higher 
than z ~ 0 have used bli nd CO surveys tuned to find CO emis¬ 
sion at z > 1. For exai^le. lAravena et al.l ( 12013) performed a Jan¬ 
sky Very Large Arra)[j blind survey towards a candidate cluster 
at z ~ 1.5, from which they obtained two detections and placed 
constraints on the number density of bright CO(l — 0) galaxies. 
IWalter et alj ^20141) carried out a blind CO survey with the Plateau 
de Bure Interferometer and placed constraints on both the CO(l—0) 
luminosity function and ph 2 in the redshift range z ~ 1 — 3. Both 
surveys indicate that ph 2 increases by an order of magnitude from 
2 ; = 0 to 2 « 2, but systematic uncertainties still dominate the 
measurements. 

There is also a wealth of literature using dust mass as a 
proxy for H 2 mass thro ugh the dust-to-gas mass ratio depen¬ 
dence on metal l icity (e.g. Boselli_eLaL I 2 OO 2 I : ISantini et alj|201 4l: 


galaxy simulations in which H 2 formation has been implemented 


ISwinbank et al.1 1201 4 iBethermin et al 


2 OI 5 I) . Using this tech¬ 


nique, and with the availability of IR photometry for large sam¬ 
ples of galaxies from Herschel, one can obtain more statistically 
meaningful samples than is possible with direct CO imaging. So 
far results from both direct CO detections and dust-derived H 2 
masses are in good agreement. For both techniques, scaling re¬ 
lations between the H 2 mass and other galaxy properties have 
been explored locally and at high redshift. It has been estab¬ 
lished that the H 2 -to-stellar mass ratio anti-correlates with stel¬ 
lar mass and that su ch a relation is presen t from 2 = 0 to at 
least 2 ~ 2.5 (e .g. Saintonge et alj|201ll: Tacconi et al.j 2013 


Santini_et^lJ_ 20W SaintongeeLal 20131 : iBothwell et^ 


2014 


Dessauges-Zavadsky et al. I I 2 OI 4 I 1 . These studies have also conclu¬ 


sively shown that there is an overall tendency for an increasing H 2 - 
to-st ellar mass ratio w i th increasing redshift at fixed stellar mass 
(e.g. lOeach et alJl201ll : ISaintonge et al.ll2013ll . This evidence sug¬ 
gests that gas in galaxies at high redshift is denser. 

The availability of high-fidelity observations of molecular gas 
in the local and distant Universe has opened up new means of 
testing galaxy formation models and simulations. This opportu¬ 
nity has been exploited using semi-analytic models which, when 
coupled with a star formation law that relates the SFR directly to 


Fu et al.ll20ld: Lasos et al. 

201 lbl:lLa 20 s et alJl2011al;lLaaos et al.l 

2014bl;|PoDDinsetalJl2014 

). However, all these models assume a 

constant H 2 to SFR conver 
servational evidence contra 

>ion efficiency, even though there is oh- 

dieting this (e.g. Saintonge et al 

I2OI3I; 


Hydrodynamic simulations are well suited to investigate the 
relation between SFR, H 2 mass, stellar mass and other galaxy prop¬ 
erties because the gas dynamics is followed self-consistently and no 
relations between the global (galaxy-wide) and local properties of 
the gas are necessarily imposed (which is the case in semi-analytic 
models, where galaxies are unresolved). Thus, one can use obser¬ 
vations of H 2 to test whether the simulations reproduce the scaling 
relations of the gas component of galaxies and isolate the physical 
drivers of the observed relations. 

Many of the scaling relations discussed above, for example 
those between the surface density of SFR, H 2 , and other properties 
of the ISM of galaxies (such as hydrostatic pressure), have been ex¬ 
plored in simulations of small portions of the ISM and in individual 


https://science.nrao.edu/facilities/via 


[Robertson & KravtsovI 1200^: 

Gnedin et alj 

2009inGlover & ClarlJ 

I 2 OI 2 I; Christensen et al, 2012 

; IWalch et al 

2014h. For example. 


the physical drivers of the relation between t he surface den¬ 
sity o f SFR and H 2 hav e been expl ored in iFeld mann et ^ 


l 201lh. iGnedin & Kravtso^ j201lh. and I Glover & Clark! (l2012h . 


IFeldmann et al.Uoilh and lonedin & Kravtsov j 201 lh show that 


the latter relation is not fundamental in the sense that large varia¬ 
tions are expected for varying physical conditions (e.g. gas metal¬ 
licity and interstellar radiation field). In addition. lOlover & ClarH 
( I2OI2I) argued against the importance of molecular cooling (e.g. 
from H 2 and CO) as the low temperatures and high-densities 
needed to form stars are reached e ven in si t uation s where the gas is 
forced to remain atomic (see also lScha^ ( l2004h ). We can imple¬ 
ment some of the relations derived from the high-resolution simula¬ 
tions above in lower-resolution, cosmological hydrodynamic sim¬ 
ulations to explore the scaling relation between the gas properties 
and other galaxy properties. 

Recently, significant advances have been achieved in cos¬ 
mological hydrodynamical simulations, which can now produce 
a galaxy population with p roperties that are broad ly consistent 
with many observat ions (see ^gelsberger^ 'ZII 2 OI 4 I for the Illus- 


with many observat ions (seel Vog elsber 
tris simulation and lSchave et alj|2015l 


for the Evolution and As¬ 
sembly of GaLaxies and their Environments, EAGLE, simula¬ 
tions). For example, EAGLE reproduces the stellar mass function 
of galaxies at z ~ 0.1, the SFR history of the Universe, and the 
properties of star-forming galaxies, in large cosmological volumes, 
(~ (100 Mpc)®), with enough particles to resolve the Jeans scale 
in th e warm ISM . _ 

iGenel et^ ( 1201 4l) presented molecular gas scaling relations 
for the Illustris simulation, and showed that it reproduces the ob¬ 
served correlation between H 2 and stellar mass reasonably well if 
the star-forming gas (which has high enough density to allow star 
formation to take place) is used as a proxy for the H 2 mass. How¬ 
ever, estimating the H 2 fraction from the neutral gas introduces 
a large systematic uncertainty because the range of temperatures, 
densities and metallicities of the star-forming gas indicates that, 
contrary to the assumption of Genel et al., it is not always H 2 - 
dominated. _ 

The OWLS simulation suite described in ISchave et alJdTOlOl) . 
a forerunner of EAGLE, explored lower-resolu t ion and smaller cos¬ 
mological volumes than EAGLE, f Duffy et ^ l l2012al) presented a 
comparison between OWLS and the H 2 mass function and scal¬ 
ing relations for galaxies with Mh 2 > 10^° M©, which was the 
mass range probed by their simulation. The authors found that the 
H 2 content of galaxies decreases with decreasing redshift, but their 
galaxies tended to be more H 2 -rich than observed and the mass 
range probed was very limited. We extend this work using more so¬ 
phisticated, theoretically motivated ways of calculating the H 2 frac¬ 
tion, and apply them to EAGLE, which has a much larger dynami¬ 
cal range than OWLS and, unlike OWLS, reproduces key properties 
of the observed galaxy population. 

Here, we present a comprehensive comparison between EA¬ 
GLE and observation s of H 2 on galaxy sc ales. Compared with, fo r 
example, the work of IPuffv et all ( 1201 2bl) and iGenel et al.l ( |2014|) . 
we extract a more realistic H 2 fraction for individual gas particles, 
based on their gas metallicity, density, pressure, SFR and Jeans 
length. The methods to compute the H 2 fraction of gas particles that 
we use here employ simple, analytic prescriptions that are based on 
more detailed studies. This kind of parametrisation avoids the use 
of expensive radiative transfer and high-resolution simulations that 
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Molecular Hydrogen in EAGLE 3 


are required to account for direct H 2 formation within a multi-phase 
ISM. However, we note that in parallel to our simple approach, the 
development of detailed chemistry netw orks and their application 
to individual galaxies is under way fc.g. Glover & JapDsenll2007l : 
iRichings et alj|201^ : iRichings et al. l2014al) . In future, the appli¬ 


cation of such networks to entire galaxies will enable us to test how 
well the phenomenological approach that we adopt here describes 
the ISM in galaxies, and to identify regimes within which our ap¬ 
proach is inaccurate. An analysis of the atomic hydrogen content 
of galaxies in EAGLE wil l be presented i n Crai n et al. (in prep.), 
Bahe et al. (submitted) and iRahmati et al.l ( l2015h . 

This paper is organised as follows. In § [^ we give a brief 
overview of the simulation and the subgrid physics included in EA¬ 
GLE. In § 0 we describe how we partition the gas into ionised, 
atomic and molecular components, the methods we use to do this 
and their main differences, with detailed equations provided in Ap¬ 
pendix!^ Comparisons with local Universe and high-redshift sur¬ 
veys are presented in § |4| and § [5] respectively. We discuss the 
results and present our conclusions in § [^ In Appendix |B] we 
present ‘weak’ and ‘strong’ convergence tests (terms introduced by 
ISchave et alj20I^ . focusing on the H 2 properties of galaxies. 


2 THE EAGLE SIMULATION 

The E AGLE simulation su itj^ (described in d etail bv ISchave et al.l 
I2OI5L hereafter S15, and ICrain et al.l I2OI5I) consists of a large 
number of cosmological simulations with different resolution, box 
sizes a nd physical models, adopting the cosmological parameters 
of the jPlanck CollaborationI l l2014ll . S15 introduced a reference 
model, within which the parameters of the subgrid models gov¬ 
erning energy feedback were calibrated to ensure a good match to 
the 2 = 0.1 galaxy stellar mass function, whilst also reproduc¬ 
ing th e observed sizes of present-day disk galaxies. iFurlong et al.l 
( I2OI4I) presented the evolution of the galaxy stellar mass function 
for the EAGLE simulations and found that the agreement extends 
to 2 « 7. 

The EAGLE suite includes simulations adopting the refer¬ 
ence model, run within cosmological volumes of (25cMpc)®, 
(50cMpc)^ and (lOOcMpc)®, where cMpc denotes comoving 
megaparsecs. In Table (Tj we summarise technical details of the 
simulations used in this work, including the number of parti¬ 
cles, volume, particle masses, and spatial resolution. The label 
of each simulation indicates box size and particle number. For 
example, L0100N1504 is a simulation of length 100 cMpc on 
a side, using 1504^ particles of dark matter and an equal num¬ 
ber of baryonic particles. Note that the spatial resolution of the 
intermediate-resolution simulations (such as L0100N1504) was 
chosen to (marginally) resolve the Jeans length in the warm ISM 
(with a temperature of « 10"* K; see SI5 for details). In Table[T] 
pkpc denotes proper kiloparsecs. 

The EAGLE simulations were performed using an extensively 
modified version of the parallel W-bo dy smoothed particle hy- 
drod ynamics (SPH ) code gadget-3 dSoringel et alj uOOa : see 
also Ispringelll200^ for the publicly available code gadget-2). 
Gadget-3 is a Lagrangian code, where a fluid is represented by 


^ See http://eagle.strw.leidenuniv.nl and 

http://icc.dur.ac.uk/Eagle/ for images, movies and data 
products. 


a discrete set of particles, from which the gravitational and hydro- 
dynamic forces are calculated. SPH properties, such as the den¬ 
sity and pressure gradients, are computed by interpolating across 
neighbouring particles. The main modifications to the standard 
Gadget-3 code for the EAGLE project include updates to the 
hydrodynamics, as described in Dalla Vecchia (in prep.; see also 
Appendix A in S15), which are collectively referred to as ‘Anar¬ 
chy’. The impact of these changes on cosmological simulations 
are discussed in Schaller et al. (in prep.). Among the main fea¬ 
tures of Ana rchy are ( i) the pressure-entropy formulation of SPH 
described in Hogkin^ 20131) , (ii) the artificial viscosity switch of 
ICullen & Dei^rJ " 201^, (iii) an a rtificial conductio n switch de¬ 
scribed in IPricd )l2008h . (iv) a C2 IWendlandl d 19951) kernel with 
58 neighbours to interpolate SPH propert ies across neighbouring 
particl es, and (v) the timestep limiter from iDurier & Dalla Vecchia 
1 I2OI2I) . which is required to model feedback events accurately. 


Another major aspect of the EAGLE project is the use of state- 
of-the-art subgrid models that capture the unresolved physics. We 
briefly discuss the subgrid physics modules adopted by EAGLE 
in § 12.11 In order to distinguish models with different parameter 
sets, a prefix is used. For example, Ref-L100N1504 corresponds 
to the reference model adopted in a simulation with the same box 
size and particle number of L100N1504. A complete description 
of the model can be found in S15 and an analysis of the impact 
of variations in the parameters of the subgrid physics on galaxy 
properties is given in icrain et alj ) l2015l) . Crain et al. also presented 
the motivation for the parameters adopted in each physics module 
in EAGLE. 


S15 introduced the concept of ‘strong’ and ‘weak’ conver¬ 
gence tests. Strong convergence refers to the case where a simu¬ 
lation is re-run with higher resolution (i.e. better mass and spatial 
resolution) adopting exactly the same subgrid physics and param¬ 
eters. Weak convergence refers to the case when a simulation is 
re-run with higher resolution but the subgrid parameters are recal¬ 
ibrated to recover, as far as possible, similar agreement with the 
adopted calibration diagnostic (in the case of EAGLE, the 2 = 0.1 
galaxy stellar mass function and disk sizes of galaxies). S15 argue 
that simulations that do not resolve the cold ISM, nor the detailed 
effects of feedback mechanisms, require calibration of the subgrid 
model for feedback, and that such simulations should be recali¬ 
brated if the resolution is changed, because in practice changing 
the resolution leads to changes in the physics of the simulation as 
better spatial resolution allows higher gas densities to be reached. 
In this spirit, S15 introduced two higher-resolution versions of EA¬ 
GLE, both in a box of (25 cMpc)® and with 2 x 752® particles, Ref- 
L025N0752 and Recal-L025N0752. These simulations have better 
spatial and mass resolution than the intermediate resolution simu¬ 
lations by factors of 2 and 8, respectively. 

Comparisons between the simulations Ref-L025N0376 and 
Ref-L025N0752 represent strong convergence tests, whilst com¬ 
parisons between the simulations Ref-L025N0376 and Recal- 
L025N0752 test weak convergence. In the simulation Recal- 
L025N0752, the values of 4 parameters have been slightly modi¬ 
fied with respect to the reference simulation. These are related to 
the efficiency of feedback associated with star formation and active 
galactic nuclei (AGN), which are adjusted to reach a similar level 
of agreement with the galaxy stellar mass function at 2 = 0.1 as 
obtained for the Ref-L025N0376 simulation (see S15 for details). 
In Appendix IB] we present strong and weak convergence tests ex¬ 
plicitly involving the H 2 abundance of galaxies. 
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Table 1. EAGLE simulations used in this paper. The columns list: (1) the name of the simulation, (2) comoving box size, (3) number of particles, (4) initial 
particle masses of gas and (5) dark matter, (6) comoving gravitational softening length, and (7) maximum proper comoving Plummer-equivalent gravitational 
softening length. Units are indicated below the name of column. The top two simulations are of intermediate resolution, while the bottom two are high- 
resolution simulations. In EAGLE we adopt (6) as the softening length al z ^ 2.8, and (7) at z < 2.8. At z = 2.8, 2.66 ckpc = O.Tpkpc and 
1.33 ckpc = 0.35 pkpc. 


(1) 

(2), 

(3) 

(4) 

(5) 

(6) 

(7) 

Name 

L 

# particles 

gas particle mass 

DM particle mass 

Softening length 

max. gravitational softening 

Units 

[cMpc] 


[M0] 

[M0] 

[ckpc] 

[pkpc] 

Ref-L025N0376 

25 

2 X 376^ 

1.81 X 10® 

9.7 X 10® 

2.66 

0.7 

Ref-L050N0752 

50 

2 X 752® 

1.81 X 10® 

9.7 X 10® 

2.66 

0.7 

Ref-L100N1504 

100 

2 X 15043 

1.81 X 10® 

9.7 X 10® 

2.66 

0.7 

Ref-L025N0752 

25 

2 X 7523 

2.26 X 10® 

1.21 X 10® 

1.33 

0.35 

Recal-L025N0752 

25 

2 X 752® 

2.26 X 10® 

1.21 X 10® 

1.33 

0.35 


2.1 Subgrid physics included in EAGLE 


• Radiative cooling and photoheating. Radiative cooling and 
heating rates are computed on an element-by-element basis for 
gas in ionisation equilibrium exposed to the Cosmic Microwave 
Background Radiation (CMB), and the model fo r the U V and X- 
ray background radiation from iHaardt & Madaul (12001). Cooling 
tables are produced with CLOUDY (version 07.02: iFerland et ^ 
Il998h for the 11 elements that dominate the cooling rate (i.e. 
H, He, C, N, O, Ne, Mg, S, Fe, Ca, Si), which are followed in¬ 
dividually. For more details see lWiersma et alj ( l2009al) and S15. 

• Star formation. Gas particles that have cooled to reach den¬ 
sities greater than are eligible for conversion to star particles, 
where is defined as: 


« in-l -3 
riH — 10 cm 


0.002 


( 1 ) 


which has a d ependence on gas metallicity, Z, as described in 
ISchavel ll2004h and SI5. Ga s particles with densities great er than 
are assigned a SFR, m* dSchave & Dalla Vecchiall2008l) : 

rh* = mg A (IM 0 pc ) " ’ (2) 

where rrtg is the mass of the gas particle, 7 = 5/3 is the ratio of 
specific heats, G is the gravitational constant, /g is the mass frac¬ 
tion in gas (which is unity for gas particles), P is the total pres¬ 
sure. A = 1.515 X 10“"^ M 0 yr“^ kpc“^ and n = 1.4 are chosen 
to rep roduce the observed Kennicutt-Schmidt relation dKennicutJ 
Il998l) . scaled to a Chabrier initial mass function (IMF). A global 
temperature floor, Teos(p), is imposed, corresponding to a poly¬ 
tropic equation of state. 


Pocp^“% (3) 

where 7 eos = 4/3. Eql^is normalised to give a temperature Teas = 
8 X 10^ K at rtH = 10~^ c m~^, which is typical of the warm ISM 
(e.g. lRichings et aT]|20I4bh . 

• Stellar evolution and enrichment. The contribution from mass 
and metal loss from Asymptotic Giant Branch (AGB) stars, massive 
stars and supernovae (b oth core collapse and type la) are tra cked 
using the yield tables of Portinari et aljdl998h , lMarigol d200ll) . and 
iThielemann et al.l (12003). The mass and metals that are lost from 
stars are added to the gas parti cles that are within the SPH kernel 
of the given star particle (see IWiersma et ai]|2009t and S 15 for 
details). The adopted stellar IMF is that of Chabrier d2003h . with 
minimum and maximum masses of 0.1 M 0 and 100 M©. 


• Stellar Feedback. The method used in EAGLE to represent 
energetic feedback associated with star f ormation (which we refer 
to as ‘ stellar feedback’) was introduced bv IPalla Vecchia & Schavel 
( l 2012 h . and consists of a stochastic selection of neighbouring gas 
particles that are heated by a temperature of 10^ ® K. A fraction of 
the energy from core-collapse supemovae is injected into the ISM 
30 Myr after the star particle forms. This fraction, /th, depends on 
the local gas metallicity and density as. 


/th 


/th,n 


/th, 


,max 


/th, 


,min 


1 + 


{ Z / ^H.b.rth 

\0.1ZqJ UH.O 


(4) 


where Z is the star particle metallicity, nn, birth is the gas den¬ 
sity of the star particle’s parent gas particle, when the star particle 
was formed. The parameters /th.min and /th.max are fixed at 0.3 
and 3, respectively, for the simulations analysed here. In the refer¬ 
ence model, nz = Tiu = 2/ln(10) and nn.o = 0.67cm“®, while 
in the recalibrated model rin is reduced to l/ln(10) and nn.o to 
0.25 cm“®. In principle, /th.max > 1 may seem unphysical. How¬ 
ever, the mechanism of stellar feedback in EAGLE includes several 
forms of feedback in the ISM: supernovae, stellar winds, radiation 
pressure, etc. Moreover, a subgrid efficiency greater than unity can 
be justified if there are numerical radiative losses. S15 quoted the 
mean and median values of /th at 2 = 0.1 which are 1.06 and 
0.7 for the Ref-L100N1504 simulations, and 1.07 and 0.93 for the 
Recal-L025N0752, respectively. 

The metallicity dependence in Eq. |4] aims to take into account 
the faster cooling that takes place in regions of high metallicity, 
while the density dependence compensates for numerical radiative 
losses that occur at high density as a result of finite resolution (see 
ICrain et ^l2015l and S15 for the motivation of this model). The 
need for this phenomenological model comes from the fact that 
stellar feedback is not simulated from first principles. 

• Black hole growth and AGN feedback. When halos become 
more massive than 10 ^° M 0 h~^, they are seeded with black holes 
of mass lO®M 0 fi“^. Subsequent gas accretion episodes make 
black holes grow at a rate that is computed following the m odi- 
fied Bondi-Hoyle accretion rate of iRosas-Guevara et al.l ( l2013l) . In 
this modification, the angular momentum of the gas can reduce the 
accretion rate compared to the standard Bondi-Hoyle rate if the tan¬ 
gential velocity of the gas is similar to, or larger than, the local 
sound speed. A viscosity parameter, Cvisc, controls the degree of 
modulation of the Bondi rate in high circulation flows, and it is a 
free parameter of the model. The value of Cvisc in the reference 
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simulations is Cvisc = Zvr, and in the recali brated, higher reso lu- 
tion simulation it is Cvisc = 27r x 10^ (see lCrain et al.ll2015l for 
analysis of the effect of Cvisc on galaxy properties). The Edding¬ 
ton limit is imposed as an upper limit to the accretion rate onto 
black holes. In addition, black holes can grow by merging with 
other black holes. 

For AGN feedback, a simila r model to the stochastic model of 
IPalla Vecchia & Schavd ( 1201 2h is applied. Particles surrounding 
the black hole are chosen randomly and heated by a temperature 
ATagn = 10 ®'® K in the reference simulation and ATagn = 
10® K in the recalibrated simulation (see Table[T]for the prefixes). 

Throughout the paper we make extensive comparisons be¬ 
tween stellar mass, SFR, H 2 mass and metallicity. Following S15, 
all these properties are measured in 3-dimensional apertures of 
30 pkpc. The effect of the aperture is minimal, even if we double 
its size to 60 pkpc differences are within 15% and much smaller for 
most galaxies. 


3 CALCULATING HI AND H 2 GAS FRACTIONS 

We describe how we calculate the neutral gas fraction of each gas 
particle in ? l3.1l and how much of that neutral gas is in molecular 
form in § 13.21 We discuss two prescriptions for calculating the H 2 
fraction and some basic predictions of each of them. For clarity, we 
introduce the number density of hydrogen, nn = nmi -h whi + 
2 nH 2 , and of neutral hydrogen, Un = whi + 2 nH 2 • We also refer 
to the fraction /h 2 = SH 2 /Sn, where Eh 2 and En are the mass 
surface densities of H 2 and neutral hydrogen, respectively. 


3.1 Calculating the neutral fraction 

Before we calculate the mass of H 2 associated with each gas par¬ 
ticle, w e need to determine th e transition from HII to neutral hy¬ 
drogen. iRahmati et alj ( 1201 3all studied the neutral gas fraction in 
cosmological simulations by coupling them to a full radiative trans¬ 
fer calculation with the aim of predicting the neutral gas column 
densities. Rahmati et al. presented fitting functions to their results, 
which enable the calculation of the neutral fraction on a particle- 
by-particle basis from the gas temperature and density, and the 
total ionisation rate (photoionisatio n plus collisional ioni sation). 
This is presented in Appendix A of IRahmati et alj ( l2013all . Here, 
we do not repeat thei r calculations but instead refer the reader to 
IRahmati et al.l j2013al) . Using these fitting functions we calculate 
77 = rin/riH. 


3.2 Calculating the H 2 gas fraction 

In this section we bri efly describe t h e pre scriptions of 
I Gnedin & Kravtsc^ ( 1201 ll) and lKjumhoi3 j2013l) to calculate 
the H 2 gas fraction for individual gas particles, which we apply to 
EAGLE. Both prescriptions give H 2 fractions that depend on the 
gas metallicity and the interstellar radiation field. 

We use the SPH kernel-smoothed metallicities rather than the 
particle metallicity with the aim of alleviating the consequences 
of the lack of me tal mixing that arises f rom the fact that metals are 
fixed to particles (IWiersma et alj2009bl) . We assume the interstellar 
radiation field to be proportional to the SFR surface density, Esfr 
(see Appendixl^for details). 


3.2.1 The prescription o ACnedin & Kravtsotl ilOlA) applied to 
EAGLE 

I Gnedin & Kravtsc^ ( 1201 ih developed a phenomenological model 
for H 2 formation and applied it to a set of zoom-in cosmological 
simulations at high resolution that also included gravity, hydro¬ 
dynamics, non-equilibrium chemistry combined with equilibrium 
cooling rates for metals, and a 3-dimensional, on the fly, treatment 
of radiative transfer, using an Adaptive Mesh Refinement (AMR) 
simulation. The highest spatial resolution achieved in these simu¬ 
lations is 260 comoving parsecs (cpc), and the gas mass of indi¬ 
vidual resolution elements ranges from 10® Mq to 10® M© . From 
the outcome of the simulations, the authors parametrised in ana¬ 
lytic formulae the fraction of H 2 -to-total neutral gas, /h 2 . as de¬ 
pendent on the dust-to-gas ratio and the interstellar radiation field. 
We assume the dust-to-gas mass ratio scales with the local metal¬ 
licity, and the radiation field scales with the local surface den¬ 
sity of star formation, which we estimate from the properties of 
gas particles assuming local hy drostatic equilibrium ( ISchavel 200 ll : 
ISchave & Dalla Vecchiall2008l) . The equations used to implement 
this prescription into EAGLE a re given in Appendix IaI We refer 
to this prescription as ‘GKl 1 ’. [Gnedin & Drainel j2014l) (GD14) 
presented an updated version of the GKll model taking into ac¬ 
count the line blending in the Lyman-Werner bands. We find that 
this model gives results that different from the GKll prescription 
by ~ 50%, in a way that the resulting H 2 masses after applying 
the GDI4 prescription are lower than when applying the GKl 1 pre¬ 
scription. This difference decreases with redshift, and even reverses 
at z > 3.5, with the H 2 masses in the GD14 model being higher 
than in the GKll model. The differences between the GDI4 and 
GKll prescriptions are of a similar magn itude as the differences 
between the GKll and iKtumholzl ( 1201 3ll prescriptions (? 13.2.4b . 
However, applying the GD14 or GKll prescriptions leads to the 
same conclusions we present in this work, and therefore we will 
hereafter only use the GKl 1 prescription. In Appendix IAll we show 
comparisons between the GKl 1 and GD14 prescriptions. 

In order to visualise how the H 2 column density is distributed 
with respect to the stars, we show in Fig. [T] five randomly chosen 
galaxies embedded in halos of different masses. The optical images 
were created using three m onochromatic rad iative transfer simula¬ 
tions with the code skirt teaes et ^ .120111) at the effective wave¬ 
lengths of the Sloan Digital Sky Survey (SDSS) u, g and r filters. 
Dust extinction was implemented using the metal distribution di¬ 
rectly from the simulations, and assuming that 30% of the metal 
mass is locked up in dust grains. Only material within a spherical 
aperture of radius of 30 pkpc was included in the radiative transfer 
calculation. The method will be presented in detail in Trayford et 
al. (in prep.). 

The variety of morphologies in the EAGLE simulations is well 
captured in Fig.[T](see also Fig. 2 in S15), where we see early-type 
galaxies (top row), SO-like galaxies (second row) and galaxy merg¬ 
ers (fourth row). In the case of the early-type galaxy example in the 
top row of Fig.[T] one can see that the H 2 is relatively concentrated 
compared to the stars, and is also connected to a dust lane. This 
is in qualitative agre ement with spatial ly resolved observations of 
early-type galaxies dYoung etalj|201 lh. in which H 2 is typically 
concentrated in the centre and in a relatively thin disk. 

3.2.2 The prescription o AKrumholA llOlA) applied to EAGLE 

lKrumhoi3 (1201 3l) developed a theoretical model for the transition 
from HI to H 2 , that depends on the total column density of neu- 
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Figure 1. Visualisation of 5 galaxies at 2 = 0 in Ref-L025N0376, which 
have been randomly chosen (and oriented) from different ranges of subhalo 
masses (indicated at the bottom of the left panels). Maps are coloured by 
H2 column density (calculated using the GKll prescription; left panels), 
and the right panels show the stellar light based on monochromatic u, g and 
r-band SDSS filter means and accounting for dust extinction (see text for 
details). Colours in the H2 column density maps are as in the colour bars at 
the top of the top, where column densities are in units of cm“^. Particles 
are smoothed by Ickpc in the Aijj, maps. Map sizes are 0.1 X 0.1 pMpc^. 


tral hydrogen, the gas metallicity and the interstellar radiation field. 
The fraction /hj is determined by the balance between photodis¬ 
sociation of H 2 molecules by the interstellar radiation field and the 
formation of molecules on the surface s of du st grains. 

A key property in the lKrumhol3 ( 1201 3l) model is the density 
of the cold neutral medium (CNM). At densities nn > 0.5 cm“®, 
the transition from HI to H 2 is mainly determined by the mini¬ 
mum density that the CNM must have to ensure pressure balance 
with the warm neutral medium (WNM, which is HI dominated). 
The assumption is that the CNM is suppor ted by turbulence, w hile 
the WNM is thermally supported (see also lWolfire et alj|2003l) . At 
riH < 0.5 cm“® the transition from HI to H 2 is mainly determined 
by the hydrostatic pressure, which has three components: the self¬ 
gravity of the WNM (oc Ehi )5 the gravity between the CNM and 
WNM (oc EhiSh 2 )> and the gravity between the WNM and the 
stellar plus dark matter component (oc EniEgd, where E^d is the 
surface density of stars plus dark matter). Note that the exact value 
of riH at which the transition between these two regimes takes place 
is a strong function of gas metallicity. The equations used to imple¬ 
ment this presc ription into EAGL E are given in Appendix We 
will refer to the lKrumholzl ( 1201 3h prescription as ‘K13’. 


3.2.3 Resolution limit on the H 2 galaxy masses 

We define the resolution limit of the intermediate-resolution simu¬ 
lations in terms of the subhalo masses rather than H 2 masses and 
present the resolution analysis in Appendix IS We find that ha¬ 
los with Mxot > 10^^° M 0 , where Mxot is the dark plus bary- 
onic subhalo mass, are converged in their H 2 abundances in the 
intermediate-resolution simulations. 


3.2.4 Differences between the GKll and K13 prescriptions 

The lower two panels of Fig.|^show the predicted ratios of the total 
H 2 mass to total neutral hydrogen gas mass (atomic plus molecular) 
as a function of stellar mass when the GKll (top panel) and K13 
(bottom panel) prescriptions are applied. There is a weak positive 
correlation between the H 2 to total neutral hydrogen mass ratio and 
stellar mass in both the Ref-L100N1504 and Recal-L025N752 sim¬ 
ulations when the K13 prescription applied, but such a trend is not 
seen in the case of the GKll prescription. The width of the distri¬ 
butions (as shown by the and ST*** percentiles) seen in both 
the Ref-L100N1504 and Recal-L025N752 simulations when the 
GKl 1 prescription is applied is much larger at Mgteiiar > M© 
than when the K13 prescription is applied. Most of the dispersion 
seen in Fig. |2] is due to gas metallicity. For example, for the sim¬ 
ulation Ref-L100N1504 with the GKll prescription applied and 
at Msteiiar ~ 10® Mq, the metallicity increases from 0.1 Zq at 
Mns/Mneutrai « 0.05 to 1Z© at Mns/Mneutrai « 0.8 (see 
also Fig. IA4I for the effect of metallicity on /h 2 )• Observations 
favou r a positive correl ation between Mh 2 /^neutral and stellar 
mass dLerov et al.ll2008l) , similar to the predictions of the Recal- 
L025N752 simulation. The very weak dependency predicted by the 
Ref-L100N1504 simulation contradicts observational results. We 
find that this is due to the gas metallicities in the Ref-L100N1504 
simulation being higher than in the Recal-L025N752 simulation 
and observations dSchave et al.l201^ . We discuss this in 1) 14.1.11 
The top panel of Fig. shows the ratio between the H 2 that is 
locked up in star-forming particles to the total H 2 , which also in¬ 
cludes the contribution from non-star-forming particles. The GKll 
prescription applied to EAGLE results in a relative contribution to 
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Figure 2. Top panel: The ratio of the H 2 mass that is locked in star-forming 
particles to total H 2 mass, as a function of stellar mass at 2 : = 0 in the 
simulation Ref-LI00N1504 after applying the GKll (solid line) and K13 
(dashed line) prescriptions. Lines and error bars indicate the median and 
16*^ and 84*^ percentiles. Middle panel. The ratio of the H 2 mass to to¬ 
tal neutral hydrogen mass, as a function of stellar mass at 2 = 0 in the 
simulation Ref-LI00N1504 after applying the GKll (solid line) and K13 
(dashed line) prescriptions. Lines and error bars are as in the top panel. Bot¬ 
tom panel: As in the middle panel but for the simulation Recal-L025N0752. 


when the GKll prescription is applied) is of a similar magnitude 
as the observed contribution from the diffuse H 2 to the total in¬ 
ferred H 2 mass. This suggests that EAGLE can shed some light 
onto the nature of such a component of the ISM of nearby galaxies. 
We investigate this in more detail in a separate paper (Lagos et al. 
in preparation). Additional differences between the GKl 1 and K13 
prescriptions, as well as comparisons with previous work can be 
found in Appendix IaI 

Fig.l^shows that the simulations Ref-L100N1504 and Recal- 
L025N0752 display significant differences both for the GKll 
and K13 prescriptions, particularly in galaxies with stellar masses 
Afsteiiar < 10^° M 0 . This is mainly driven by the gas metallicity 
in the Recal-L025N0752 simulation being lower on average than 
in the Ref-L100N1504 simulation (see Fig. 13 in S15 for the mass- 
metallicity relations predicted by the simulations compared here). 
For example, at a stellar mass of 10® M©, Recal-L025N0752 sim¬ 
ulation predicts a gas metallicity 0.4 dex lower than the metallicity 
of galaxies of the same mass in the Ref-L100N1504 simulation, on 
average. However, for Mgteiiar > 10^® Mq differences are small. 
In ? l4.Lll we discuss the effect the gas metallicities have on the H 2 
masses of galaxies in EAGLE. 


4 ABUNDANCE AND SCALING RELATIONS OF H 2 IN 
THE LOCAL UNIVERSE 

In this Section we compare the results of EAGLE with observations 
of H 2 in the local Universe. We focus on the H 2 mass function in 
ij 14. H and on H 2 scaling relations in 8) 14.21 

Note that in previous papers using EAGLE, simulation results 
and observations are compared at 2 = 0.1 in order to best cor¬ 
respond to the median redshift of optical surveys. However, here 
we compare the EAGLE predictions with observations at 2 = 0 
because observations of CO, which here we convert to H 2 , in the 
local Universe correspond to galaxies at 2 < 0.1. 


the total H 2 mass from the H 2 locked up in non-star-forming parti¬ 
cles that is on average 20 — 40 per cent and does not show a strong 
dependence on stellar mass. However, when the K13 prescription 
is applied, this contribution decreases to 5 — 20 per cent (see the 
discussion in Appendix lAl for the contribution to the total H 2 mass 
of gas particles that are on the equation of state). This difference is 
significant, given that the errors on the median assessed via boost¬ 
rapping are very small (< 0.05 dex). The source of the difference 
between the GKl 1 and K13 prescriptions, is that the K13 prescrip¬ 
tion leads to /hj = 0 at gas surface densities Eh < 1 M© pc“® 
in ISM conditions that are MW-like, while leading to /hj = 0 at 
higher Eh when the ISRF increases. On the contrary, the GKll 
prescription leads to small, but nonzero H2 fractions at these low 
gas surface densities. 

It has been observed in nearby galaxies that there is a compo¬ 
nent of the molecular gas that is not associated with star formation 
and that is chara cterised by a high velocity dispersion con sistent 
with the HI disk jPetv et alj|200 . ICaldu-FTimo et alj|2013h . This 
component contributes « 30 — 50% of the total H 2 mass inferred 
for those galaxies. This observed diffuse component could be as¬ 
sociated with the H 2 locked up in non-star-forming particles in 
EAGLE, as their contribution to the total H 2 in galaxies (at least 


4.1 The H 2 mass function 

The curves in Fig. [3 show the 2 = 0 H 2 mass function calculated 
using either the GKll or K13 prescriptions, and the mass func¬ 
tion of the hydrogen component of the star-forming gas, for the 
simulation Ref-L100N1504. These are compared to observations 
at 2 « 0, which are shown as symbols with \a errorbars. We can 
see that the star-forming gas mass is a poor proxy for the H 2 mass. 
The median ratio between H 2 mass and the hydrogen component 
of star-forming gas is ~ 0.7, but fluctuations around this value are 
large, with the lO*** and 84*'’ percentiles being 0.33 and 1.25, re¬ 
spectively. We qualitatively discussed the large systematics intro¬ 
duced by the assumption that the star-forming gas is a good proxy 
for H 2 mass in § 1 in the context of previous work (e.g. lGenel et al.l 
l2014h . and here we have quantified them. 

As we discussed in § 1, H 2 is not directly observed, but con¬ 
straints on the H 2 mass function at 2 = 0 have been provided by 
surveys of CO(l-O) in local galaxies. The CO(l- 0 )-to-H 2 conver¬ 
sion factor is the largest uncertainty in the observations. The con¬ 
version factor, X, is defined as 


-= X X 10 

cm~^ 


( -^CO(l-O) \ 

VKkms-i ) ’ 


(5) 


where Xh 2 is the column density of H 2 and Jco is the integrated 
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Figure 3. The mass functions of H 2 (thick lines) and of the hydrogen com¬ 
ponent of the star-forming gas (thin line) in the simulation Ref-L100N1504 
at z = 0. For H 2 we show the mass functions after applying the GKll 
(black line) and K13 (red line) prescriptions. The part of the mass functions 
that is affected by low number statistics (less than 10 objects per dex bin) is 
shown as dashed lines, while the part that is affected by resolution is shown 
as dotted lines (see Appendix [^ for details). Observational estimates cor¬ 
respond to the iKeres et al.l (200^ .B-band (triangles) and 60 /im (circles) 
selected samples of galaxies, adopting the Milky-Way like H2-to-CO(l-0) 
conversion factor X = 2 as universal (s ee Eq.[5]for the definition of X). W e 
also show the observational estimates of lObreschkow & Rawlingg ^20091) in 
square symb ols, in which a lum inosity-dependent X was applied to the ob¬ 
servations of iKeres et aljj2003) . The star-forming gas mass function differs 
significantly from the H 2 mass functions obtained using either the GKl 1 or 
the K13 prescriptions at > 5 x 10® Mq, and can therefore not be 

used as a reliable tracer of the H 2 mass function. The H 2 mass function 
predicted by EAGLE is in good agreement with the observations. 


CO(l — 0) line intensity per unit surface area. Systematic varia¬ 
tions of X are both theoretically predicted and observationally in¬ 
ferred. For example, variations are expected as a function of gas 
metallicity and the interstellar radiation field, such that X increases 
in low metalli city and/or intense interstellar radiation field envi- 
ronrn e nts (e.g. Bell et alJ [2007 : Meiierink et alj 200^: Bavet et al 


2 OO 9 I: IShettv etalJl201lt 


s aravanan 


et alJ I 2 OI ?: iBisbas et al 


20151) . In addition to this, different values of X have been measured 
using different methods. For example, in the Milky Way, variations 
in the metho d for measuring X lead to values that differ by a factor 
of « 3 Isee lSolatto et al.ll2013l for a recent review). This system¬ 
atic error clearly dominates the uncertainty in the H 2 mass over the 
uncertainties in the measurement of the CO emission line. We keep 
these systematics in mind when comparing the EAGLE predictions 
with the observations. 

Eor the entire range of observed H 2 masses, 10^ M© < 
Mh 2 ^ 10^° M 0 , the H 2 mass function of EAGLE, calculated 
either using the GKll or the K13 prescriptions, is in good agree¬ 
ment with the observations. We s ee the agreement is the best with 
the H 2 mass function deriv ed by Obreschkow & Rawlingsl l l2009l) 
from the CO(l-O) surveys of iKeres et al. j2003h using a luminositv- 


dependent H2-to-CO(l-0) conversion factor. Compared to the H 2 
mass functions derived using a fixed X — 2, we find some ten¬ 
sion at Mhj ^ 5 X 10® Mq. However, the differences are well 
within the systematic uncertainties of the observational measure¬ 
ments, which are well illus trated by the difference betwee n the 
observational infe rences of lObreschkow & Rawlingg (I2009I) and 
iKeres et"^ ( l2003h . 


4.1.1 Variations due to gas metallicity 


S15 compared the mass-metallicity relation (MZR) predicted by 
EAGLE with observations and showed that the intermediate- 
resolution simulation Ref-L100N1504 exhibits ISM metallicities 
that are higher than the observations by up to 0.3 dex at a stel¬ 
lar mass of < 10® Mq. The higher resolution simulation, Recal- 
L025N0752, on the other hand, displays ISM metallicities that are 
lower than the lower-resolution simulations at fixed stellar mass, 
and are in good agreement with the observations. 

The gas metallicity is important for two reasons. Eirst, the 
GKll and K13 prescriptions have an explicit dependence on gas 
metallicity, as dust is the catalyst for H 2 formation, and the abun¬ 
dance of dust is assumed to be proportional to the gas metallicity. 
Second, the formation of CO, and therefore the reliab ility of CO as 
a trac er of H 2 , depends on the gas metallicity (e.e. iBolatto et al.l 
l2013h . Observations sug gest a strong depen dence of X on gas 
metallicity. Eor example, Boselli et alj ( ItOOtI) reported the follow¬ 
ing relation between X and Z, 

logio(2f) = 0.5t°:'2 - 1 . 02 tH 5 logio(Z/Zo). (6) 


In order to quantify the effect that the predicted higher metal¬ 
licities have on the H 2 content of galaxies in EAGLE, we first re¬ 
peat the calculation presented in ? 13.21 but forcing galaxies to fol¬ 
low the observed MZR. In practice, we compare the neutral gas 
mass-weighted metallicity of indiv idual EAGLE galax ies with the 
observed median gas metallicity of IZahid et alJ l l2014l) at the same 
stellar mass and estimate how much higher the predicted metallic¬ 
ity is compared to the median in the observations. If the difference 
is greater than 50% of the observed median, we rescale the metal¬ 
licities of all gas particles by the factor needed to bring the galaxy 
metallicity into agreement with the observations, and recalculate 
the H 2 fraction and the 30 pkpc aperture H 2 mass. The 50% fac¬ 
tor is included to allow for a dispersion around the median of the 
same magnitude as the la distributions reported in the observa¬ 
tions. We convert the metallicities of Zahid et al. to our adopted 
metal abundances: 12 + logio(0/H)Q = 8.69 and Zq = 0.0127. 
Eor stellar masses Mgteiiar < 10® M© we extrapolate the MZR 
of Zahid et al. to get an ‘observed’ metallicity at that stellar mass. 
The H 2 mass functions using the GKl 1 prescription before and af¬ 
ter the rescaling of gas metallicities are shown for the simulations 
Ref-LIOON 1504 and Recal-L025N0752 in the top panel of Eig.|4] 
The differences in the H 2 mass functions before and after 
the metallicity rescaling are mild. The metallicity rescaling has 
no visible effect on the H 2 mass function of Recal-L025N0752, 
while for Ref-LIOON 1504, the number density of galaxies with 
Mh 2 ^ 5 X 10® Mq decreases slightly after rescaling. When in¬ 
tegrating over all H 2 masses, the average density of H 2 in Ref- 
L100N1504 decreases by 5% after the gas metallicity rescaling. 
This difference is well within the error bars of the observations (see 
S I5.4t . However, larger differences are obtained for the relation be¬ 
tween the Mna/ATneutrai ratio and stellar mass (Eig.O, in which 
the gas metallicity rescaling leads to Mna/Afneutrai ratios lower 
by up to Ri 0.5 dex at 10® Mq < Mateiiar 5 X 10® Mq, for both 
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Figure 4. Top panel: The H 2 mass function at 2 : = 0 for the simulations 
Ref-L100N1504 and Recal-L025N0752 using the GKll prescription, be¬ 
fore and after the gas meta llicities of galaxies in EAGLE are rescaled to 
satisfy the observed MZR of lZahid et al.H2014h . We refer to the simulations 
after rescaling as ‘FMZR’. As shown in Fig. [3] the effects of resolution in 
the Ref-L100N1504 start to appear at H 2 masses ^ 5 x 10^ Mq. The ef¬ 
fect of sampling, on the other hand, appears at number densities of —2.5 
in the Recal-L025N0752 and —4.3 in the Ref-L100N1504 (in the units of 
the y-axis). Observations are as in Fig. im The rescaling of gas metallici¬ 
ties has a significant effect at < 10® Mh 2 - Bottom panel: CO(l-O) 

luminosity function for the same models of the top panel, and applying the 
conversion factor of Eq. ito convert the predicted H 2 masses into CO(l-O) 
luminosities. Ob servations correspo nd to the S-band and 60 microns se¬ 
lected samples of iKeres et alj j2003h . as labelled. 


prescriptions. A steeper relation is then obtained, which removes 
the tension with the observations described in § 13.2.41 

As a second experiment, we use Eq. to convert the pre¬ 
dicted H 2 masses in EAGLE to CO(l-O) luminosities and con¬ 
struct the CO(l-O) luminosity function, in order to show the effect 
of a metallicity-dependent X. This is shown in the b ottom panel 
of Fig. |4l along with the observations of CO(l-O) of iKeres et al.l 


(l2003h . The level of agreement between the Ref-L100N1504 and 
Recal-L025N0752 simulations and the observations is similarly 
good to that for the H 2 mass functions, and the effect of rescaling 
the gas metallicity is also similar. This shows that the main effect of 
gas metallicities is in the H 2 formation rather than its later conver¬ 
sion to CO(l-O). Note that the CO(l-O) luminosity functions pre¬ 
dicted by the simulations Ref-L100N1504 and Recal-L025N0752 
are closer to each other than the predicted H 2 mass functions. 
This is because the simulation Recal-L025N0752 has, on average, 
lower gas metallicities than the Ref-L100N1504 simulation, which 
in practice leads to higher X values. 

The effects of the higher gas metallicities on the H 2 masses of 
galaxies is mild and therefore we continue to use the predicted gas 
metallicities (rather than the rescaled ones) for the rest of the paper. 


4.2 H 2 Scaling relations 


Measurements of the H 2 and HI mass content, as well as other 
galaxy properties, are available for relatively large samples of lo¬ 
cal galaxies (running into the hundreds), enabling the character¬ 
isation of scaling rel ations between the cold gas and the stel- 
lar mass conten t (e.g . [Bothwell et al.ll2009l : ICatinella et al.ll2oT^ : 
ISaintonge et al.l 1201 ll : iBoselli et alJ 2014h . Recently, H 2 and HI 
have been s tudied in homogeneou s sa mples of relatively ma ssive 
galaxies by ISaintonge et ah (1201 ih an d ICatinella et alj ( l2010t) . re¬ 
spectively, with the aim of measuring the fundamental relations be¬ 
tween the stellar content of gala xies and their cold gas mass. The 
survey of ISaintonge et alj 1 I 2 OI ll) . the CO Legacy Database for the 
GALEX Arecibo SDSS Survey (COLD GASS) is the first stel¬ 
lar mass-limited CO survey. The strategy used was to select all 
galaxies with Msteiiar > 10^° Mq at z < 0.05 from the SDSS 
and follow-up a subsample of those at millimetre (mm) wave¬ 
lengths to detect CO(l-O). Saintonge et al. integrated sufficiently 
long to enable H 2 gas fractions, Mns/ATateiiar > 0.015, at stel¬ 
lar masses Msteiiar > 10 ^°'® Mq, to be detected, or H 2 masses 
Mh 2 > 10®'® M© in galaxies with stellar masses, lO*^® Mq < 
Mateiiar < 10^®'® M© . The simple and well-defined selection func¬ 
tion of COLD GASS makes this survey an instructive testbed for 
EAGLE. 


Fig.H] shows the H 2 fraction, Mh 2 /(ATh 2 + ATsteiiar), as a 
function of stellar mass for galaxies with Msteiiar > lO*^® M© 
and H 2 gas fraction s above the sensitivity limit defined by 
ISaintonge etZIdTOllh for COLD GASS at z = 0 in EAGLE. Note 
that both observations and simulation results account for molecu¬ 
lar hydrogen only (without Helium correction). We show the me¬ 
dian and 16*® and 84*® percentiles for galaxies in COLD GASS, 
for which the CO(l-O) emission line is detected as squares with 
error bars, while we show the upper limits for the non-detections 
as arrows. Note that uncertainties in the measurement of the inte¬ 
grated flux in the CO(l-O) emission line are small (approximately 
1%). Here the systematic errors are dominant, and they are mainly 
related to how the CO flux is converted into the H 2 mass, as we 
discussed in § 14.11 

From Fig. we see that EAGLE predicts that the ratio 
Mh2/(Mh 2 + Mateiiar) decreases with Msteiur- The amplitude of 
this relation is slightly higher, by < 0.1 dex, when the GKll pre¬ 
scription is applied than when the K13 prescription is applied. The 
high-resolution simulation Recal-L025N0752 predicts a trend that 
is similar to that in Ref-L100N1504, although with systematically 
higher Mh 2 /(Mhs + Msteiiar) (except for the largest galaxies, for 
which we have only a small number in the Recal-L025N0752 sim¬ 
ulation). This is due to the slightly higher SFRs and H 2 masses 
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Figure 5. Molecular hydrogen fraction, + Mgteilar). as a 

function of stellar mass for the simulations Ref-L100N1504 (solid line and 
filled region) and Recal-L025N0752 (dashed line with error bars) for galax¬ 
ies with Mgteilar > 10 ^*^ Mq and with H 2 masses above the sensitivity 
limit of COLD GASS aX z = 0. Here we use the GKll prescription to 
calculate H 2 masses. We also show the simulation Ref-L100N1504 when 
the K13 prescription is applied (long-dashed line) for the same redshift and 
stellar mass selection. Lines show the medians of the simulated galaxies 
for which the H 2 mass exceeds the COLD GASS detection limit, while the 
filled region shows the 16**' and 84*** percentiles for the Ref-L100N1504 
simulation when the GKl 1 prescription is applied. The dispersions in the 
Ref-L100N1504 simulation with the K13 prescription and in the Recal- 
L025N0752 simulation are of similar size. Squares with error bars show 
the median and 16*** and 84**^ percentiles of galaxies in the COLD GASS 
survey that have the CO(l-O) emi ssion line detected, w hile an'ows show 
upper limits for the non-detections iSaintonge et al.l201 ih . The predictions 
from EAGLE are in very good agreement with the data. 


predicted in the higher resolution simulation at a fixed halo mass. 
S15 showed that the Ref-L100N1504 simulation predicts SFRs that 
are lower than the observations by ~ 0.1 — 0.15 dex for galaxies 
with Msteiiar > 10*° Mq, while the Recal-L025N0752 simulation 
predicts higher SFRs, in better agreement with the observations. 
Such differences also appear in Fig. |5] where the median in the 
Ref-L100N1504 simulation is ~ 0.1 dex lower than the median H2 
fraction in COLD GASS. This offset is smaller than the scatter in 
the observations. Thus, the agreement seen here is fully consistent 
with the analysis of the specific SFRs in S15. The differences be¬ 
tween the Ref-L100N1504 and Recal-L025N0752 simulations are 
also similar to those discussed for the H 2 -subhalo mass relation in 
Appendix IB] 

The predictions from the EAGLE simulations shown in Fig.|5] 
agree very well with the observations. Galaxies in COLD GASS 
for which the CO(l-O) emission line is detected lie on the main 
sequence of star-forming galaxies (in the SER-stellar mass plane), 
which is also the case for galaxies in EAGLE that have H 2 masses 
above the detection limit of COLD GASS. This is a significant out¬ 
come, since none of these observations were used to constrain the 
calibrated parameters governing the efficiency of star formation and 
feedback, nor those governing the partitioning of gas into ionised, 
atomic and molecular components in EAGLE. This is the first time 



Figure 6. The SFR as a function of the molecular hydrogen mass at z = 0 
for galaxies with ATstellar > 10 *** Mq and with H 2 gas fractions above the 
sensitivity limit of COLD GASS in the simulations Ref-L100N1504 (solid 
line and filled region) and Recal-L025N0752 (dashed line) with the GKll 
prescription, and the simulation Ref-L100N1504 with the K13 prescription 
(long dashed line), where lines show the medians. Observations are as in 
Fig.E] The simulations and observations are in good agreement. 


that a hydrodynamic simulation performs so well when confronted 
by observational inferences of H 2 . 

In Fig. [6] we show the SFR as a function of the H 2 mass for 
EAGLE galaxies with Msteiiar > 10*° Mq and H 2 gas fractions 
above the sensitivity limit of COLD GASS at 2 = 0 and for the 
COLD GASS observations. In EAGLE there is a tight correlation 
between these two quantities, because star formation and most of 
the H 2 coexist in the same gas phase (see Fig. This is also 
the reason why the relation between SFR and Mh 2 is the tight¬ 
est of all the scaling relations discussed in this section. Simula¬ 
tion Recal-L025N0752 displays a lower normalisation than Ref- 
L100N1504 (with the same prescription applied) at H 2 masses 
Mh 2 ^ 10°'® Mq by a factor of ~ 1.5, consistent with the typ¬ 
ically greater H 2 mass associated with galaxies at a fixed stel¬ 
lar mass in the Recal-L025N0752 simulation, as discussed in Ap- 
pendix|2 Note that the tight relation between the SFR and the H 2 
mass is not surprising, because the star formation law adopted in 
EAGLE relates the SER of particles with their gas mass, above a 
given density threshold. However, the component of the H 2 mass 
that is non star-forming gas, in addition to the dependence on gas 
metallicity, result in differences between the different simulations 
with the K13 or GKll prescriptions to look different in the SFR- 
Mh 2 relation (with differences being of ~ 0.2 dex), even though 
the star formation law adopted is the same in every case. 

The relation between SER and the H 2 mass predicted by EA¬ 
GLE agrees with the reported detections in COLD GASS. Galaxies 
for which the CO(l-O) emission line is not detected have signifi¬ 
cantly lower SERs. EAGLE predicts that the H 2 masses of these 
galaxies should be an order of magnitude (or more) below the de¬ 
tection limit of COLD GASS. This is because EAGLE predicts that 
the tight relation between SFR and H 2 extends down to very low 
values in both quantities. This can be seen in Fig.|3 where we lower 
the threshold in Mh 2 /Msteiiar to Mh 2 /Msteiiar > 0 (while main¬ 
taining the threshold in stellar mass). 
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Figure 7. As Fig. |6] but for galaxies with /Msteiiai > 0 and 
^stellar > 10^° Mq. Here the filled region con'esponds to the 2.5**' and 
97.5**' percentiles for the simulation Ref-L100N1504 with the GKIl pre¬ 
scription applied. To the observations shown in Fig.[6l we add those of the 
ATLAS®*^ survey iYoung et al.l201 ll i lPavis et alj201^ . as labelled, which 
sample a similar range of stellar mass as COLD GASS, but go much deeper 
in the CO(l-O) emission line. The main difference is that ATLAS®*** only 
includes early-type galaxies, while COLD GASS does not distinguish be¬ 
tween different morphological types. Thus, ATLAS®** cannot be directly 
compared to EAGLE (given the selection) but one can conclude that about 
17% of the early-type galaxies in ATLAS®** are > 2a from the median 
relation predicted by EAGLE. 


In Fig.[3we also show the observations of the ATLAS®** sur¬ 
vey lICaDDellari et al.ll201 ll) . which covers a similar range in stellar 
mass as COLD GASS but only includes early-type galaxies (el¬ 
lipticals and lenticulars). The complementary aspect of ATLAS®** 
relative to C OLD GASS is tha t it detects much fain ter CO(l- 
0) fluxes (see I Young et al.ll201ll and lYoung et alJl2oI3) . SFRs in 
ATLAS ®** come from a c ombination of UV and mid-IR photome¬ 
try (see lDavis et aLll^oTil for details). We find that Ri 17% of the 
early-type galaxies in ATLAS®** lie > 2a from the median relation 
in EAGLE, pointing to a possible conflict between EAGLE and 
the observations. However, this cannot be confirmed here because 
we would need to make a thorough comparison with ATLAS®** 
by morphologically selecting galaxies to be early-type in EAGLE. 
This is beyond the scope of this paper. If the discrepancy was real, 
it could be due to the star form ation law used in EAGLE, which is 
constructed from the observed iKennicutl (Il998h relation between 
the surface densities of SFR and gas. Galaxies in ATLAS®** as well 
as the u pper limits from C OLD GASS (shown in Fig. [7} are offset 
from the lKennicud ( Il998h relation in the sense that their efficiency 
of star formation is inferred to be lower (see also conclusions of 
iDavis et^l2014h . 

Fig. [H shows the Hg depletion timescale, tb _2 = MH 2 /SFR, 
as a function of the stellar surface density at 2 ; = 0 for the simu¬ 
lation Ref-L100N1504 with the GKl 1 prescription and the median 
and 16*** and 84*** percentiles of the observations of COLD GASS. 
The H 2 depletion timescale is defined as the time that would be 
required to exhaust the current H 2 mass if the SFR were constant. 
The stellar surface density is defined in EAGLE and in the obser¬ 
vations as fii, = (Msteiiar/2)/(7rr5o), where Msteiiar and rso in 


Figure 8. The molecular hydrogen depletion timescale, = 

Mjjj/SFR, as a function of the stellar' mass surface density, /r* = 
Afgtellar /at z = 0 (see text for details on how /r* is calculated in ob¬ 
servations and EAGLE), for galaxies in EAGLE with ATgtellar > 10 *® Mq 

and H2 gas fractions above the sensitivity limit defined for COLD GASS). 
Simulations and observations shown are as in Fig.|^ 


EAGLE are the stellar mass in the subhalo and the projected radius 
that contains half of that stellar mass, respectively. In COLD GASS 
the stellar mass is derived from the broad-band photometry and rso 
is the radius encompassing half of the Petros ian flux in the 2 -band, 
which is measured in circular apertures (see lSaintonge et ^|2012| 
for details). We show the median and Ifi*** and 84**' percentiles of 
galaxies in COLD GASS for which the CO(l-O) emission line is 
detected as squares with error bars, while upper limits are shown as 


Even though the EAGLE predictions lie in a similar region 
to the observations in Eig.[8] there are some interesting discrepan¬ 
cies. It has been suggested from the observations shown in Fig. [8] 
that there is a tendency for tho to increase with increasi ng /r* 
dSaintonge et alJ[20l^ : lMartig et al.|[20r3l : IPavis et alj|20i^ . This 
trend is apparent if upper limits are considered, but less clear if 
only detections in COLD GASS are considered. It has been sug¬ 
gested that such a trend can arise if the central stellar surface den¬ 
sity increases the surfac e density threshold for gas instability, there- 
fore loweri ng the SFR dSaintonge et alj|2012 : iMartig et al.ll2013l : 
iDavis et aDl2014l) . in other words, as dynamically driven star for¬ 
mation suppression. Such a trend is not seen in the reference sim¬ 
ulation, Ref-L100N1504, when the GKll prescription is applied, 
but instead the predicted relation is consistent with a constant ths ■ 

To gain more insight into the relation between ths and /r*, we 
again relax the limit imposed on Mh 2 /ATsteiiar to Mh 2 /ATsteiiar > 
0 , while maintaining the cut in stellar mass, Msteiiar > 10*® M©. 
The results are shown in Fig. where we also include the observa¬ 
tions of the ATLAS®** survey, in which /r* is measured from the K- 
band derived stellar masses and the half-light radius in the r-band. 
Since most early-type galaxies in ATLAS®** do not have signifi¬ 
cant ongoing star formation, measuring rso in the r-band is similar 
to measuring it in the infrared (as was done for COLD GASS). 

The observations of ATLAS®** confirm the existence of a pop¬ 
ulation of galaxies with larger values of rH 2 than the galaxies in 
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Figure 9. As Fig. [8]but for galaxies with Mhj/ j\<fstellar > 0. Here the 
filled region con'esponds to the 2.5*^ and 97.5**' percentiles of the simula¬ 
tion Ref-L100N1504 with the GKll prescription applied. As in Fig.|2]we 
add the observations of ATLAS®***, for which /r* is calculated using the r- 
band half-mass radius. ATLAS®*** cannot be directly compared to EAGLE 
(given the selection) but one can conclude that about 7% of the early-type 
galaxies in ATLAS®*** are > 2a from the median relation predicted by EA¬ 
GLE. 


COLD GASS with CO(l-O) line detections. The trends predicted 
by EAGLE are as in Fig.jS] although the dispersion around the me¬ 
dians becomes much larger at /i* > 10* ® Mekpc”^. This means 
that some galaxies in EAGLE that have large values of /r* also 
have high values for thj , consistent with the observations, but not 
the majority of them. However, we find that only 7% of early-type 
galaxies in ATLAS*® are > 2a from the median relation predicted 
by EAGLE, and again, we cannot confirm this is a real issue in the 
simulation given that we are not selecting galaxies in EAGLE to be 
early type. Although we have not applied the same definition of /r* 
used in the observational data, we find that the trends of Fig.[8]are 
not strongly sensitive to the definition of /r*. For example, similar 
trends and values for /i* and thj are found if we take the half-mass 
radius measured in apertures of SOpkpc and 20pkpc. 

We conclude that EAGLE captures the main scaling relations 
between stellar mass, SFR, stellar mass surface density and H 2 
mass of galaxies that are on the main sequence of star formation. 
Some discrepancies are present when we consider galaxies with 
high values of Ta 2 , which in EAGLE are rare but observationally 
represent a large fraction of the early-type galaxy population in the 
local Universe. These tests are important since none of these ob¬ 
servations were used to fix the free parameters in EAGLE. A study 
of the radial profiles of these quantities is ongoing and will be pre¬ 
sented in a future paper (Lagos et al., in preparation). 


5 THE EVOLUTION OF THE H 2 ABUNDANCE OF 
GALAXIES 

In this section we present the analysis of the evolution of the H 2 
mass function (1) 15. It , H 2 scaling relations (l; l5.2t . the stellar masses 
and SFRs of galaxies selected by their H 2 mass (§ ES, and the 
global H 2 density (? l5.4b . over the redshift range 0 2 : 5. 


5.1 The evolution of the H 2 mass function 


At _2 = 0 t he H 2 mass function has be en measured bv iKeres et al.l 
( I 2 OO 3 I) and lObreschkow et alj l l2009ah . At z — 1.5 and 2 = 2.7 
constraints on the H 2 mass function ha ve been provided fro m a 
blind molecular em ission line survey by IPecarli et alj (12014l) and 
IWalter et al.l (|2014|) . Although the uncertainties in these measure¬ 
ments are significant (see Fig. [JoJ, these are currently the best 
available measurements of the H 2 mass function at these redshifts. 
Other constraints on the H 2 mass functions at 2 « 1 — 2 have 
been obtained by following-up galaxies at mm wavelengths that 
have spectroscopic o r photometric redshifts. This is the case for 
lAravena et al.l ( l2012h . who performed a JVLA survey towards a 
candidate cluster at 2 ~ 1.5. Aravena et al. obtained two detections 
and were able to place constraints on the number density of bright 
CO(l — 0) gala xies. Another inter esting measurement at 2 « 1.5 
was provided bv IPaddi et al.l ( l201(]h . who detected CO(2— 1) emis¬ 
sion lines for six galaxies from a sample of colour-selected star¬ 
forming galaxies. 


Fig. [To] shows the H 2 mass function of the simulation Ref- 
L100N1504 using both the GKll and K13 prescriptions for sev¬ 
eral redshift bins ranging from 2 = 0 to 3. The H 2 mass function 
shows an increase in the number density of galaxies of all masses 
from 2 = 0 to 2 ~ 1.2 — 1.5. At 2 > 1.5 the trend reverses and 
the number density of galaxies starts to decrease with redshift. The 
driver of this reversal is that at 2 > 1 the gas metallicity decreases 
faster than at 2 < 1, which hampers H 2 formation. In addition, the 
median interstellar radiation field in galaxies increases monotoni- 
cally with redshift (at least up to 2 = 4.5), hampering H 2 formation 
even further. We come back to this point in 1) 15.41 Interestingly, the 
differences at 2 = 0 between the GKll and K13 prescriptions ap¬ 
plied to EAGLE disappear at 2 « 1 — 2, and at 2 « 2 — 3 the 
K13 prescription gives slightly higher abundances of galaxies with 
Mh 2 < 10® Mq than the GKll prescription. This is because the 
GKll prescription gives lower H 2 fractions than the K13 prescrip¬ 
tion in the case of intense interstellar radiation fields (Go > 100 ). 
The latter values of Gq are typical at 2 > 1.5 (see bottom panel 
of Fig. [T^ and Appendix for an analysis of the differences be¬ 
tween the H 2 masses predicted when applying the GKll and K13 
prescriptions). 

We find that EAGLE predictions at 2 « 1.2 (middle panel 
of Fig. [To), with either the GKll and K13 prescripti ons applied, 
are w it hin the uncertainties of the observations o f IPaddi et aP 
ll 2010 l) . lAravena et all jiol^ and IWalter et all ( |20l4 . At z m 2.5 
(right-hand panel of Fig. llOt . the predicted abundance of galax- 
ies with Muo > 1 0 ^° Mq is low compared to the constraints of 
IWalter et al.l ( 20141) . particularly if we compare with the observa¬ 
tional constraints at 2 = 3. However, the redshift range of that data 
point in the right panel of Fig. [T^ is large, 2 « 2 — 3.1, which can 
partially explain the offset with the predictions. Another important 
u ncertainty is the CO -H 2 conversion factor: X — 2 was assumed 
in IWalter et alj ( l201^ and lPecarli et al.l ( l2014h . If this conversation 
factor were on average lower in galaxies with large abundances o f 
H 2 , as suggested by th e the oretical studies of Lagos et ^ ( l2012l) . 
iNaravanan et alj (12012l) a ndlPopping et al.l (|2014|) . the H 2 masses 
would be overestimated by [Walter et al. 1 2014h . improving the con¬ 
sistency with the EAGLE results. 

Note that the small number of galaxies with Mua ^ lO*^® M© 
is not an effect of the limited box size, as we find from compar¬ 
ing the Ref-L100N1504 and Ref-L050N0752 simulations which 
produce a very similar number density of galaxies with Mb _2 ^ 
10*® Mq and a very similar value of the cosmic mean H 2 density. 
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Figure 10. The H 2 mass functions for the reference simulation Ref-L100N1504 using the GKll (solid lines) and K13 (dashed lines) prescriptions applied, 
at z = 0 and 2 : = 0.5 in the left-hand panel; z = 1.2 and z = 1.9 in the middle panel; and 2 = 2.3 and 2 = 3 in the right-hand panel. The pai1 affected 
by resolution is approximately 5 x 10^ Mq (see Appendix [B^. Observatio ns shown in the left pan el are as in Fig. [3 Th e solid grey regions in the 

middle and right-hand panels con'espond to the blind CO measurements presented of lWalter et alJ i2Q14l) an d lPecarli et alJi2Ql3) . at 2 1.5 and 2 ^ 2.7, 

respectively. The upper side of the grey regions show the number density obtained if all the high-quality line candidates are considered real, while the bottom 
side of the boxes show the number density if only secure detections are cons idered. The error b ars i n the boxes indicate t he Poisson errors on the number 
densities. In the middle panel we also show the observational inferences from lPaddi et alj j2Q10t) and lAravena et al.l i2Q12h . as labelled, which ai'e based on 
CO emission line follow-up of 2 1.5 galaxies or fields with spectroscopic or photometric redshifts. 


(Fewer than 10 objects per dex bin in Fig.llOlcorresponds to a num¬ 
ber density of 10~^'^ in the units of the y-axis.) 

The EAGLE results in Fig. \M are slightly different from 
those given by semi-analytic models. The latter predict the peak 
of the numb er density of galaxies with Mho 10^° M q to be at 
2 ~ 2 — 2.5 dObreschkow et al.l2009bl : lLagos et al.l2012h . A possi¬ 
ble explanation for this is that the models that have been preferred 
in semi-analytic work to calculate t he ratio of H 2 -to-HI gas m ass 
(for example the empirical model of iBlitz & RosolowskvIbOO^ do 
not consider gas metallicity. Future observations can distinguish 
between these different predictions, for example with telescopes 
such as the Atacama Large Millimeter Array (ALMA). 


5.2 H 2 scaling relations at high redshift 

Fig.im shows the ratio -|- Msteiiar) and the FI 2 de¬ 

pletion timescale, ths = Mna/SFR, as a function of redshift, 
calculated for galaxies on the main sequence of star formation 
and with stellar masses Msteiiar > 5 x 10® M© in the Ref- 
L100N1504 simulation using the GKl 1 and K13 prescriptions. We 
define the main sequence of galaxies as ±0.5 dex around the me¬ 
dian logic (SFR) at a given stellar mass, c alculated from the simu ¬ 
lation outpu t s. Obs er vational data are from Saintonge_et^ 
lOeach et al.l d201 1). Magdis et alj ( 1201 2h . Saintonge et al. 
and Tacconi et alj 2013 1. and also refer to quantities measured for 
main sequence galaxies over a range of stellar masses similar to the 
one chosen here. The definition of the main sequence of galaxies 
in the observational data is also based on the median SFR of galax¬ 
ies of a given stellar mass. However, since observations are flux 
limited, the possibility exists that the observational data are biased 
towards higher SFRs. 


EAGLE predicts the ratio Mh 2 /(Mh 2 ± Msteiiar) to increase 
by 0.5 — 0.7 dex from z = 0 to 1 and to evolve very little from 
2 « 1 to 4. These trends are also observed. The GKl 1 prescription 
applied to the Ref-L100N1504 simulation predicts a slightly lower 
median value of Mhj /(Mhj ±Msteiiar) at 2 > 2.5 than is the case 
for the K13 prescription. The reason is that the GKl 1 prescription 
is more sensitive to intense interstellar radiation fields, Gq > 100, 
than the K13 prescription, and the latter values are typical for z >2 
(see bottom panel of Eig. m. 

There is a slight discrepancy of a factor of « 1.5 between the 
predicted ratio Mh 2 /(Mh 2 ± Msteiiar) and the observationally in¬ 
ferred ones at z « 2 — 3, regardless of whether the GKll or the 
K13 prescription is applied. It is unclear at this stage whether this 
offset is a problem for EAGLE, since the observations have so far 
been biased towards gas-rich systems (as they only quote detec¬ 
tions, except for the values in the local Universe given by COLD 
GASS). However, the offset between the predicted and observed ra¬ 
tio Mh 2 /(Mh 2 ± Msteiiar) is, in any case, within the scatter seen 
in observations. 

The increase of the H 2 mass fraction with redshift is driven by 
the greater fraction of gas particles that reach the regions of high 
H 2 fractions (T < 10"^ ® K and P/ks ^ 5 x 10® cm~®). Such a 
trend can be seen from the evolution of the Mh 2 -weighted pressure 
in the top panel of Fig. d The trend shows that at higher redshift 
more H 2 is locked up in particles of higher pressure than at 2 = 0. 

The bottom panel of Fig.ll llshows the evolution of the deple¬ 
tion timescale, rH 2 , measured for the same main sequence galaxies 
used to calculate the evolution of Mh 2 /(Mh 2 ± Msteiiar) shown 
in the top panel. We find that rH 2 decreases with redshift by a fac¬ 
tor of ~ 10 between z = 0 and z = 4, with the GKll prescription 
predicting a slightly larger decrease than the K13 prescription. Such 
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Figure 11. Top panel: The mean H 2 gas fraction, Mh 2 / (-^H 2 +-^stellar) 
as a function of redshift for the Ref-L100N1504 simulation using the GKl 1 
(solid line) and the K13 (dashed line) prescriptions. The mean H 2 gas frac¬ 
tion was calculated in EAGLE using only main sequence galaxies with 
Mgteiiar > 5 X 10^ Mq. The observational data also correspond to the me¬ 
dian H 2 gas fraction for main sequence galaxies at different redshifts. Data 


was taken fromISaintonae et alj i201lb,lGeach et al 

j201lh.lMa2dis et alj 

^20121). ISaintonae et alj J2013h and iTacconi et alj 

201 3|). Magdis et al. 


presented a detection and a non-detection of the CO(3-2) emission line 
of galaxies at 2 : 3, and therefore the combined estimate presented 

here is only an upper limit. Bottom panel: The H 2 depletion timescale, 
MH 2 /SFR, as a function of redshift. The obseiwational data are as in the 
top panel except for the measure ment at 2 : 2.4 which is a combination of 

the lTacconi et alj i2013h and the lSaintonge et al] j2013h measurements, as 
was presented by Saintonge et al. 
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Figure 12. Top panel: The median of the Mh 2 -weighted pressure, P, when 
the GKll (solid line) and the K13 (long dashed line) prescriptions are 
applied to the simulation Ref-L100N1504, as a function of redshift. The 
SFR-weighted P (dotted line) is also shown as a function of redshift in 
the same simulation. Pressure is expressed in terms of P fcg where /cb is 
Boltzmann’s constant. Bottom panel: The median Ssfb/Ssfr,0 ’ where 
SsFR,o = 10~^ M(7) yr~^ kpc~^ is the SFR surface density in the solar 


neighborhood iBonatto & Bic j|201 ll) . calculated using all the star-forming 
particles at each time, as a function of redshift. 


evolution indicates that EAGLE galaxies at higher redshift are more 
efficient at forming stars. The non-monotonic behaviour of the H 2 
fraction with redshift, seen at the top panel of Fig. m affects TH 2 
by accelerating the decrease with redshift at 2 : > 2.3. This can be 
seen from fitting a power-law relation before and after 2 = 2.3: at 
2 ^ 2.3, th2 oc (1 -I- 2)“'^'^, while at 2 > 2.3, rH2 oc (1 + z)~^. 

To help visualize the reasons behind the increase in the effi¬ 
ciency of star formation with redshift, we show in the top panel 
of Fig. [m the SFR-weighted and Mh 2 -weighted pressures for 
the Ref-L100N1504 simulation, which we refer to as (P)sfr 
and (P)mh 2 respectively. Applying the GKll or K13 prescrip¬ 
tions in this simulation does not make any significant difference to 
('P)mh 2 . Another important observation is that both (P)sfr and 
(T’)mh 2 increase with redshift, but (P)sfr does it to a greater ex¬ 
tent. Since the star formation law adopted in EAGLE has a power- 
law, n = 1.4 (see Eq|2j, forming stars at higher pressure implies 
having higher gas surface density and therefore higher SsFu/Sgas. 
Although (P)mh 2 increases, its weaker evolution compared to 
(P) sFR helps th, to d ecrease faster with redshift. 

iLillv et akl ( 1201 3h have suggested that the observations shown 
in Eig. [TT] are consistent with the “equilibrium model”, where the 
evolution in the specific SER and the H 2 depletion timescale are 
mainly driven by the evolution of the gas reservoir (neutral gas in 
the galaxy). We find that besides the gas reservoir, an additional 
parameter plays an important role: the typical pressures at which 
star formation takes place and at which Ft 2 forms. In other words, 
galaxies move along the Kennicutt-Schmidt relation towards higher 
gas surface densities as redshift increases. The change in the FI 2 
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depletion timescale is therefore a result of adopting a superlinear 
star formation law (i.e. power-law index > 1 in Eq.|3 and the tran¬ 
sition from HI- to H 2 -dominated gas evolving towards higher mean 
ISM pressures as the redshift increases. Note that our conclusion is 
subject to the assumption that the Kennicutt-Schmidt relation holds 
at redshifts higher than 0. This is a reasonable assumption since 
observations so far show that ga laxies at z > 0 fol l ow the same 
Kenn icutt-Schmidt relation (e.g. lOenzel et al.ll201(il : iGenzel et al.l 

I 2 OI 3 I) . 

The trends obtained in EAGLE for the evolution of H 2 gas 
fractions and H 2 depletion timescales are qualitatively the same 
as those inferred from observations. However, observational re- 
sult s have so far presen ted contradictory conclusions. Eor exam¬ 
ple, iGenzeT^t^yoT^ presented CO and dust mass observations 
for ~ 500 galaxies from 2 « 0 to 2 « 3, and concluded that 
the H 2 depletion timescale is only a weak function of redshift, 
thj oc (1 -I- 2 )“°'®, while EAGLE favours a stronger dependence, 
tb _2 oc (1 -f for the Ref-L100N1504 with the K13 prescrip¬ 
tion applied and ths oc (1 -I- with the GKl l prescription 

applie d. However, recent ALMA observations by IScovilk et al.l 
l l2015h favour a much stronger redshift dependence of thjEI with 
z ~ 1 galaxies having ths lower by an order of magnitude than 
2 = 0 galaxies, in g o od agr eement with the EAGLE results. More¬ 
over, l^oville et alJ yOl^ found no differe nce in th^ between 
galaxies on and off the main sequence, while iGenzel et al.l ( 1201 Sl) 
found an order of magnitude difference. This argues for further ob¬ 
servations and the need for galaxy samples with simpler selection 
functions that can then be applied to simulations to make a one-to- 
one comparison (as we did in S I4.2l for COLD GASS). 


5.3 The properties of H 2 mass-selected samples 

An interesting question is what are the stellar masses and SERs 
of galaxies selected by their H 2 mass? This question is relevant 
for planned blind CO surveys as an important aspect of these sur- 
veys is finding opt ical and near-IR counterparts (see for example 
IPecarli et'^l2014l) . We show in Fig. [T3] the galaxy stellar mass 
functions for all galaxies (dashed lines) and for the subsample that 
have Mb_2 > 10^ Mq (solid lines) at 0 < 2 < 3. Typically, galax¬ 
ies with large reservoirs of H 2 also have high stellar masses. How¬ 
ever, at 2 < 0.5 we observe an increasing number of massive galax¬ 
ies with Mh 2 <10® M 0 , which is clear from the difference in nor¬ 
malisation between the galaxy stellar mass function of all galaxies 
and of those with Mh 2 > 10® Mq in the top two panels of Fig.fTS] 
At 2 > 1, galaxies with Mgteiiar > 10^® Mq almost always have 
Mhj > 10® Mq. For a comparison between the pr edicted galaxy 
stellar mass function in EAGLE and observations see lFurlong et al.l 
( |2014|) . 

Fig. [M] shows the SFR function of all galaxies (dashed lines) 
and of those with Mh 2 > 10® Mq (solid lines). The H2 mass- 
selected sample is a very good tracer of the most actively star¬ 
forming galaxies at any time. For example, at 2 = 0 galaxies with 
Mh 2 > 10® Mq have SFR> 1 Mq yr“^. The average SFR of the 
H 2 mass-selected sample increases with redshift, not only because 
galaxies have on average higher SERs at fixed stellar mass, but also 
because the ratio between the H 2 mass and the SFR decreases with 


® IScoville et alj l2015h measure ISM masses rather than H 2 masses. How¬ 
ever, from the very high H 2 gas fractions measured in high -2 galaxies, one 
can safely assume that the ISM masses of high -2 galaxies are dominated by 
H 2 , on average. 



9.0 9.5 10.0 10.5 11.0 11.5 

l09lo(^stella/^s) 


Figure 13. Galaxy stellar mass functions of all galaxies (dashed lines) 
and those having Mh 2 > 10® Mq (solid lines) at different redshifts, 
as labelled in each panel, for the simulation Ref-L100N1504 using the 
GKll prescription. The H 2 mass-selected sample traces galaxies with stel¬ 
lar masses, Mgteilar <) 10^® Mq, but at 2 < 0.5 about half of those 
massive galaxies have Mh 2 < 10® Mq. 


redshift (see bottom panel of Fig. [ID- This means that the H 2 mass 
at a fixed SFR decreases with redshift. The outcome of this is that 
the typical SERs of galaxies with Mh 2 > 10® Mq at 2 = 3 are 
SFR> 10 Mq yr“^ 


5.4 The evolution of Gh 2 

Fig.ES] shows the evolution of f 2 H 2 = PH 2 /Pcrit, where pcrtt 
is the critical density of the Universe at the given redshift, for 
Ref-L100N1504 using both the GKll and K13 prescriptions. The 
global density of H 2 was calculated using the H 2 abundance in 
galaxies measured in a 30 pkpc aperture. If we instead use all the 
gas particles that are bound to the sub-halo of a galaxy, we find a 
value higher by 1.5 — 3%, which shows that most of the H 2 in the 
universe resides close to galaxy centres. We also show for refer¬ 
ence the evolution of the SFR cosmic density (dotted line) in the 
Ref-L100N1504 simulation. 

The GKll prescription applied to EAGLE predicts a peak at 
a slightly lower redshift, 2 « 1.2, compared with 2 « 1.4 for 
the K13 prescription. The latter results in an increase of a factor 
of « 4.3 between 2 = 0 and 2 ~ 1.4, while GKll yields an 
increase of a factor of « 3.5 over the same redshift range. The 
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Figure 14. SFR functions of all galaxies (dashed lines) and those having 
^b .2 > ^0 (solid lines) at different redshifts, as labelled in each 

panel, for the simulation Ref-L100N1504 using the GKl 1 prescription. The 
H 2 mass-selected sample traces galaxies with SFR> 1 Mq yr~^ at 2: = 0, 
while at 2 = 3 the typical galaxies present in H 2 mass-selected samples 
have SFR> 10 Mq yr“^ at 2 : = 3. 


difference is driven by the GKl 1 prescription being more sensitive 
to intense interstellar radiation fields, as discussed in g 15:21 The 
peak of f 2 H 2 is located at a redshift which in both cases is lower 
than the peak of th e SFR density in EA GLE, at z « 2 (see dotted 
line in Fig.llSIand lFurlong et alj[2014l for a detailed analysis on 
the SFR density evolution). The reason for this offset is the faster 
decrease in gas metallicity at z > 1 and the monotonic increase of 
the median Esfr with redshift (see the bottom panel of Fig. [T^ . At 
z < 1, the FI 2 mass-weighted metallicity decreases as logio(2) oc 
(1 -F 2 )“° ®, while at z > 1 the decrease accelerates and scales as 
logio(2) oc (1 -F 2 )“^'®. Note that we still see a rapid increase 
of f 2 H 2 between 0 and 1 even though some of the conditions for 
FI 2 formation become less favourable (lower gas metallicities and 
higher interstellar radiation fields). This is due to an increase in 
the pressure of the ISM in galaxies and to an increase in the total 
neutral gas content. 

We also s how i n Fig. JBL the observatio na l est imates 
of Keres et al] (l2003h . lObreschkow & Rawlings! ( l2009h and 
I Walter et all ( 1201 4 The error bars on the Walter et al. inferences 
are as described in §111] Compared to current observational esti¬ 
mates, EAGLE captures the increase in the H 2 content of galaxies, 
but predicts a lower r2H2 in the redshift range 2 = 1 — 3, although 



Figure 15. ’ defined as the ratio of the H 2 density to the critical density 

of the universe, r 2 H 2 = PH 2 /Pcrit ^ a function of redshift for the simula¬ 
tion Ref-L100N1504 using the GKl 1 (solid line) and the K13 (dashed line) 
prescriptions. For reference, we show the SFR cosmic density evolution as 
dotted line, offset by 1.7 dex so that it lies in a similar range than • 
Obser vations at 2 : = 0 bv iKeres et al.l i2003l) and lObreschkow & Rawlin^ 
J2009tl are shown as a filled circle and star, respecti vely, with fq- error b ars, 
while the observations from the blind CO survey of lWalter et al.l i2014ll are 
shown as solid and hatched regions. The boxes show the contribution of 
the CO blind detections without extrapolating to account for the undetected 
sources below the detection limits in the millimeter (which would translate 
into H 2 masses < 5 X 10® Mq). The en'or bars in the boxes indicate Itr 
Poisson errors. 


within the error bars of the observations. Some possibilities to ex¬ 
plain this tension are: (i) EAGLE is truly predicting a lower value 
of flH 2 than in the real Universe, and (ii) the observational bias 
nH 2 towards higher values. Concerning (i), this could be related to 
the SER cosmic density in EAGLE be ing below the ob servations by 
a facto r of Ri 2 jFurlong et al.ll2014l) . Related to (ii), I Walter et^ 
(I2OI4I) adopted a value of X to convert their CO measurements to 
H 2 similar to the Milky-Way value. However, theoretical evidence 
points to a lower value of X bei ng more reasonable for these high- 
redshift galaxies. Eor example, iNaravanan et al.l ( 12012h show that 
the value of X decreases with increasing velocity width of the CO 
line. Since the CO velocity widths of these high-redshift galaxies 
are much larger than for 2 = 0 galaxies, one would expect lower 
values of X. These systematics in the observations need to be care¬ 
fully addressed before we can conclude that the tension between 
EAGLE and the observational inferences of Oh 2 is a true short- 
com ing of the simulatio n. 

[Walter eTaP ( I20l4 argue that the tension seen between sev¬ 
eral predictions of Oh 2 in the literature and their observational es¬ 
timates are, in reality, worse than it seems from comparisons such 
as the one presented here. This is because their estimate of nH 2 ac¬ 
counts only for detections (secure and candidates) and it does not 
include the contribution from galaxies with Mh 2 ^ 5 x 10® Mq. 
However, to assess whether this really worsens the tension with 
theoretical predictions, one needs to estimate how much of flH 2 
is contributed by galaxies with different values of Mh 2 or other 
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Figure 16. r2H2 as a function of redshift for the simulation Ref-L100N1504 
using the GKll prescription separated into the contribution from galaxies 
with different H 2 masses (top panel), star formation rates (middle panel) 
and stellar masses (bottom panel), as labelled. The solid line in each panel 
corresponds to the total H 2 density also shown as a solid line in Fig. [B] 
Obseiwations are as in Fig. |15l and are shown in all panels. 


properties. We show in Fig. M the value of flH 2 for subsamples 
of galaxies in EAGLE selected by their H 2 mass, SER and stellar 
mass. Galaxies with Mh 2 > 10^ Mq dominate the global flH 2 at 
2 < 2.5, which means that we expect the H 2 mass function inte¬ 
grated over the massive end only (above the break) to recover most 
of the H 2 in the universe for these redshifts. 

In terms of stellar mass and SER, flH 2 is dominated by galax¬ 
ies with: 

• SER > 10 Mq yr-i at z > 0.1, 

• Afstellar > 5 X lO'^ Mq at Z < 2 

• Afstellar < 5 X lO'^ Mq at Z > 2. 

The dominant contribution of galaxies with Mgteiiar < 5 x 10® Mq 
at 2 > 2 is unlikely to be a consequence of the limited simulated 
volume (100 cMpc box size). In fact, flH 2 convergences to better 
than 15% in the simulations Ref-L100N1504 and Ref-L050N0752. 

These results are sligh tly different to thos e found in semi- 
analytic models. For examnle fLagos et aPd^Mah find that galaxies 
with 1 Mq yr“^ < SFR < 10 Mq yr“^ make a larger contribu¬ 
tion to flH 2 than those with SFR > 10 Mq yr“^. These differences 
are interesting and will be tested by future observations. 

The galaxies that dominate flH 2 in EAGLE are readily ob¬ 


servable in current optical and near-IR surveys (for example the 
CANDELS survey: Icrogin et a0201 ih . This suggests that a cheap 
and effective observing strategy to unveil f 2 H 2 is to follow up ex¬ 
isting optical and near-infrared surveys to measure CO, rather than 
blind CO surveys that are more expensive and harder to charac¬ 
terise. In follow-up strategies, molecular emission line surveys can 
be guided by well-known UV, optical, IR and/or radio surveys, 
with well identified positions and spectroscopic or photometric red- 
shifts. Even in the case that only galaxies with SFR> 10 Mq yr“^ 
are followed up, we expect them to uncover 50% or more of flH 2 . 
placing strong constraints on theoretical models of galaxy forma¬ 
tion. 


6 CONCLUSIONS 

We have presented a comprehensive comparison between two hy¬ 
drodynamic simulations of the EAGLE suite and surveys of molec¬ 
ular gas in galaxies at low and high redshift. To determine H 2 
masses in simulated galaxies, we first extracted the fraction of neu¬ 
tral gas (atomic and molecular) in each gas particle using the fit¬ 
ting functions of iRahmati et al.l ( 12013ah . which were derived ap¬ 
plying radiative transfer algorithms to cosmological simulations in 
order to calculate local neutral fractions as a function of the colli- 
sional and photoionisation rates, temperature and gas density. Af¬ 
ter estimating the fraction of the gas that is neutral, we proceeded 
to calculate the fraction of that gas that is molecular. We applied 
two schemes to EAGLE to calc ulate the H 2 fraction of ind ividual 
gas particl e s, tho se advanced bv iGnedin & Kravtsov! ( 1201 ih and by 
lKrumhoi3 (l2013h . The first one, GKll, provides fitting functions 
for the H 2 fraction that depend on the dust-to-gas mass ratio and 
the local interstellar radiation field. These fitting functions were 
calculated from a set of small, high-resolution cosmological simu¬ 
lations, which included on-the-fly radiative transfer, a simple model 
for H 2 formation and dissociation, where the main catalyst for H 2 
formation are dust grains. The second prescription, K13, is a theo¬ 
retical model where the H 2 fraction is calculated assuming pressure 
balance between the warm and cold neutral media, and depends 
again on the dust-to-gas mass ratio and the local interstellar radia¬ 
tion field. After computing the H 2 fraction in each gas particle, we 
used apertures of 30 pkpc to calculate the H 2 mass in individual 
galaxies. 

We compared the predicted H 2 masses with observations, and 
our main conclusions are: 

• The H 2 galaxy mass functions at 2 = 0 predicted by EA¬ 
GLE for the two prescriptions tested here are very similar to 
ea ch other and i n good agre ement with the observational esti mates 
of iKeres et ^ (l2003b and lObreschkow & Rawlingsl ( l2009lt (see 
Fig.lDl. We showed that the total star-forming gas mass in EAGLE 
becomes a poor proxy for the H 2 mass for Mb _2 it 5x10® Mq . We 
find that H 2 masses become affected by resolution below Mb _2 ~ 
5 X 10^ Mq, and for Mstellar < 10® Mq. 

• W e made an extensive comparison with the COLD GASS 
survey dSaintonge et al.ll201 ih . applying the same selection crite¬ 
ria as were applied to the observations: Mstellar > 10^® Mq and 
Mh 2 /Mstellar > 0.015 in galaxies with Mstellar > 10 ^° ® Mq or 
H2 masses > 10®'® Mq in galaxies with 10^® Mq < Mstellar < 
IOIO .6 Mq 'We found generally good agreement in the scaling 
between the H 2 mass and other properties such as stellar mass 
(Fig.S, SFR (Fig. 1^ and stellar surface density (Fig. [ 8 ]l. An in¬ 
teresting discrepancy was seen in the scaling between the SFR and 
the H 2 mass. Observations reveal a population of galaxies that have 
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low SFRs, but a considerable amou nt of H 2 (for example the early- 
type galaxies in IPavis et al. I I 2 OI 4 see Fig. |7]l. These galaxies are 
offset from the main relation traced by spiral galaxies. In EAGLE 
we find that all galaxies follow the latter relation regardless of 
their morphology, with little scatter, and therefore no galaxies lie 
in the region of low SFR for a fixed H 2 mass. This is because at 
small scal e s the s tar formation law in EAGLE follows the standard 
iKennicutj ( 1 19981) star formation relation, which has a higher nor- 
malisatio n than the star fo rmation relation inferred for early-type 
galaxies jPavis et al.ll2014l) . 

• We presented the SERs and stellar masses of H 2 mass-selected 
samples of galaxies in EAGLE. The most H 2 -rich galaxies are also 
the most actively star-forming galaxies at any redshift (Fig. d, 
while also being the most massive at 2 ; > 0.5 (Fig. ll3t . At 2 < 0.5 
a significant fraction of galaxies with Mgteiiar > Mq are H 2 - 
poor (Mh 2 < 10® Mq). 

• EAGLE predicts (for both the GKll or K13 prescriptions ap¬ 
plied) that Gh 2 increases by a factor of ~ 4 from 2 = 0 to 
2 ~ 1.3 — 1.5, followed by a decrease towards higher redshifts 
(Fig.d . In EAGLE the peak redshift of flH 2 is lower than the peak 
redshift of the SFR density, 2 « 2. We find that this is due to the gas 
metallicity decreasing faster at 2 > 1 compared to the decrease at 
2 < 1, and the SFR surface density increasing with increasing red¬ 
shift, on average. Both trends hamper the formation of H 2 at 2 > 1. 
We find that the EAGLE H 2 mass functions and flH 2 agree well 
with the observational inferences in the redshift range 0 < 2 < 1.5, 
but underestimate the H 2 abundance at 2 < 2 < 3 (see Fig.fTOt. 
However, it is unclear how much of the latter discrepancy is driven 
by systematic errors in the observations. We also find that Gh 2 in 
EAGLE is dominated by galaxies with Mh 2 > 10® Mq at 2 < 2.5, 
SFR> 10 Mq yr“^ at 0.1 < 2 < 5 and Msteiiar > 5 x 10® M© 
at 2 < 2 (Fig. [16}. The latter implies that an efficient strategy for 
unveiling f2H2 is to use current optical and near-IR surveys with 
spectroscopic or good photometric redshifts and accurate galaxy 
positions and to follow them up using millimeter telescopes target¬ 
ing the CO emission. 

The agreement between the properties of molecular gas in the 
EAGLE simulations and observations at low and high redshift is 
remarkable, particularly when one considers that all the free pa¬ 
rameters of the subgrid physics models in the simulation were fixed 
by requiring a match to the observed galaxy stellar mass function 
and sizes at the present day. No parameters were adjusted to match 
any of the H 2 observations discussed in this paper. We now plan 
to study the distribution of H 2 within galaxies in EAGLE and to 
compare with resolved observations, in an effort to further our un¬ 
derstanding of how the physical processes modelled in EAGLE de¬ 
termine the properties of dense gas in galaxies. 
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APPENDIX A: THE TRANSITION FROM NEUTRAL TO 
MOLECULAR HYDROGEN 

In this Appendix we provide the equations used to apply the GKl 1 
and K13 prescriptions to EAGLE. 


Al The prescription of I Gnedin & Kravts^ j201lh applied to 
EAGLE 


In this section w e de scribe the analytic formulae that 
iGnedin & Kravtso^ ( 1201 ih obtained for computing the H 2 
fraction and the way we apply them to EAGLE. The H 2 fraction is 
calculated as 



where, Eh 2 is the surface density of H 2 and En is the surface den¬ 
sity of neutral hydrogen. The parameter Ec is defined as 

2 1 

Ec = 20 Mo pc-= -- . , (A2) 


where Dmw is the dust-to-gas mass ratio in units of the Milky 
Way dust-to-gas mass ratio, which we take to be Dmw = ZjZ^ 
under the assumption that the dust-to-gas ratio scales linearly 
with gas metallicity, Z, and that this relation is universal, Zo = 
0.0127 is the solar metallicity dAllende Prieto et al] 1200 ih . and 
Go is the interstellar radiation field (defined in the wavelength 
range 912A to 2400A) in units of the Habing radiation field, 
1.6 X 10“® ergcm“^ dHabindl 19681) . Radiation at these fre¬ 
quencies is produced mainly by OB stars. In EAGLE we adopt 
a universal IME and therefore the rate at which OB stars are 
formed is proportional to the SER, such that we can express Gq as 
Gq = Esfr/Esfr,o, where Esfr.o = 10“® M p^ yr~^ kpc~^ is 
the SF R surface density in the solar neighborhood dBonatto & Bical 
l201lh . Einally, A is defined as 


A = ln 


1 


^gD 


3/7 

MW 



(A3) 


Here, the function g is 

1 -I- as -I- 

g = -r--, 

1 -I- s 

where a and s are defined as 


s 


Q 


0.04 

Dt + Dmw ’ 


5 go/2 

1 + (Go/2)2 • 


Here D* is defined as 


(A4) 


(A5) 

(A6) 


D. = 1.5 X lO’^ln (^1-b [3Go]^ '^) . (A7) 

Eollowing ISchavel d200lh . we estimate the local neutral hy¬ 
drogen surface density, En, by multiplying the volume density of 
hydrogen by the neutral fraction and the Jeans length, Aj: 


E. 


gpH Aj = gpn 


Cs 

x/CpH 


(A8) 


Here, Cs is the effective sound speed w hich is a function of the 
press ure, P, and the gas density, p (see ISchave & Dalla Vecchi3 
|2008| for details). Erom EQS. IAll and lA8l we can estimate Eh 2 . 

We apply the above prescriptions to individual gas particles, 
where Esfr = Psfr Aj, ps fr is the SER density (compute d from 
the pressure using Eq. 12 of ISchave & Dalla Vecchiall2008h . Aj is 
calculated as in Eq IASI Un is the particle density of neutral hydro¬ 
gen and Z the particle smoothed metallicity. 
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For non star-forming particles (those that have densities 
below the density threshold for star formation, see Eq. [T](, 
pgPB, = 0, so we ado pt the UV photoionisation background of 
lldaardt & Madaul 1 2001). which is a function of redshift, to deter¬ 
mine Go ■ Rahmati et alj (120130) showed that star-forming galaxies 
produce a galactic scale photoionisation rate of « 1.3 x 10“^® s“^, 
which is of a similar magnitude than the UV background at 2 = 0, 
and smaller than it at 2 > 0. This favours our approximation, which 
uses the UV photoionisation background as the relevant photo¬ 
ionisation rate for non star-forming particles. 

Note that, since Egs.lATI to lATI are fits to the simulation set 
of iGnedin & Kravtsov! ( 1201 ll) . there are no free parameters we can 
adjust. Instead, we apply the prescription as presented in Gnedin 
& Kravtsov directly to EAGLE. We also imple mented the updated 
versio n of the GKll prescription, described in iGnedin & Drain3 
(I2OI4I) and found that H 2 masses were typically lower when the 
updated prescription is applied (which we refer to as GD14) than 
with the GKl 1 prescription. The top panel of Fig. lAl [ shows the H 2 
masses predicted by the GD14 prescription as a function of the H 2 
masses resulting from applying the GKll prescription, at differ¬ 
ent redshifts. The differences seen in the H 2 masses with the GKl 1 
and GD14 prescriptions depends on redshift, because the correction 
introduced by GD14 to the GKll prescriptions changes the depen¬ 
dency of /h 2 on the ISRF. In the bottom panel of Fig. lAll we show 
the effect the changed on Mhs have on Gh 2 • Although there is a 
change in overall normalisation, the behaviour of the two models 
is quite similar: there is a peak in Gh 2 at 2 ~ 1-2 followed by a 
decline at higher redshift. 


A2 The prescription of lKrumholzl ^20131) applied to EAGLE 

iKrumhoij ( 1201 3h considered an ISM composed of a warm neutral 
medium (WNM) and a cold neutral medium (CNM) which are in 
pressure equilibri um. This impl i es a m inimum density of the CNM, 
first described bv lWolfire et alJ (l2003h . which is 


^CNM,min ~ 31 Go 


Z^/Z 


i + 3.i(G;2d/Ct)“-3®5 


(A9) 


where Gq is the interstellar radiation field in units of the Habing 
radiation field, Zd and Z are the dust and gas phase metallici- 
ties, respectively, and (;jt is the ionisation rate due to cosmic rays 
and X-rays. Assuming that the co smic rays density and the UV in¬ 
tensity scale with the local SFR, lKrumhoi3 ( l2013l) approximates 
Go/Ct = 1- Similarly, the metallicity of the dust and gas phase are 
assumed to be equal, Z^jZ = 1. In 2-phase equilibrium the CNM 
can have de nsities between n cNM.min and « 5ncNM,min. With 
this in mind iKrumholzl ( 1201 3h writes a fiducial CNM density for 
2-phase equilibrium 


OQ /^1+3-1Gmw®A „„-3 /Aim 

ncNM, 2 p « 23 Go I -- I cm . (AlO) 

Note that we slightly changed the notation of lKrumholzl ( 1201 3h to 
make it consistent with the notation used in s IaTI 

In the regime where Go 0, ncNM, 2 p in Eg. lAlOl also tends 
to 0, which would imply that the pressure o f the CNM tends to 0 
too. This regime is unphysical, and therefore iKrumhoQ ( l2013l) ar¬ 
gues that the thermal pressure of the gas becomes relevant in this 
regime and sets a minimum C NM density to main tain hydrostatic 
balance, ncNM,hydro (see also lOstriker et alJl^l^ . The latter de¬ 
pends on t he maximum tempe rature at which the CNM can exist 
(« 243 K: IWolfire et'^l20C)3l) . 7cNM,max, thc matter density (of 




Figure Al. Top panel: The H 2 masses in the Ref-L100N1504 simulation 
when the GD14 prescription is applied compared to when the GKl 1 is ap¬ 
plied, at different redshifts, as labelled. The diagonal solid line shows the re¬ 
lation -ATh 2 (GD14) = Mh 2 (GKll). Lines with errorbars show the me¬ 
dian and 16^*^ and 84^^ percentiles. Bottom panel: The evolution of flH 2 
in the Ref-L100N1504 simulation with the GKll and GD14 prescriptions 
applied, as labelled. Symbols are as in Fig. [B] 


stars and dark matter), psd, and the neutral ga s surfa ce density, En 
(which we compute using EQ. IA8t . lKrumhol3 (1201 3l) writes: 


tt-CNM,hydro 


Pth 

1.1 Tcnm ,max 


(All) 


where 


Pth 


ttGE^ [ / 32Cda/wCipsd 

+ V -^GEl- 


(A12) 


Here a « 5 represents how much of the midplane pressure sup¬ 
port comes from turbulence, ma gnetic fields and cosm ic rays, com¬ 
pared to the thermal pressure dOstriker et alj|201^ . (^d ~ 0.33 
is a numerical factor that depends on the shape of the gas sur- 
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face isodensity contour, /w = 0.5 is the ratio between the mass- 
weighted mean square thermal velocity dispersion and the square 
of the sound speed of the warm gas (th e value adopted here orig¬ 
inally comes from lOstriker et alJ[201011 and Cw = 8kms~^ is 
the s ounds speed in the warm neutral medium (e.g. iLerov et^ 
l2008t) . The value of the gas density in the CNM is then taken to 
be n CNM = maxfncN M. 2d, rtCNM. hydro). 

iKrumhofj ( 1201 3h defines a dimensionless radiation field pa¬ 
rameter: 


X 



n-CNM 
10 cm~3 


) 


-1 


(A13) 


and writes the ratio between the H 2 density and the total neutral 
hydrogen density as 


/h2 


1 - 0.75 s/(l-f 0.25 s), s<2 

0, s ^ 2 


where 


(A14) 



ln(l-F 0.6 x-b 0.01 x") 

0.6 Tc 

(A15) 

0.066/cDmw (m^pc-s)' 

(A16) 


(A17) 


Here /c is a clumping factor representing the ratio between the sur¬ 
face density of a gas complex and the surface density of the dif¬ 
fuse mediu m. This value is sugg e sted t o be /c « 5 on scales of 
~ 1 kpc bv lKrumholz & McKee! ll200.'5ll based on observations of 
nearby galaxies. Since the spatial resolution of EAGLE is of the 
order of 1 kpc at low redshift (see Table[TJ, we adopt /c = 5. 

As in S IAll we apply Egs. lAlOl to lAlTI to individual gas parti¬ 
cles to obtain the H 2 mass per particle, and then calculate the total 
H 2 mass within 30 pkpc of the centre of individual subhalos. 

A3 Differences between the GKll and K13 prescriptions 

Fig. IA2I shows the ratio between the H 2 mass calculated using 
the GKll prescription and the mass calculated using the K13 
prescription, Mh 2 (GK11)/Mh 2 (K13), as a function of the neu¬ 
tral gas metallicity at z = 0 in the simulation Ref-L100N1504. 
The largest differences between Mh 2(GK11) and Mh 2(K13) 
are obtained in the low-metallicity regime, Z < 0.5 Zq, where 
Mh 2 (GK11)/Mh 2(K13) > 5. However, the variation around 
the median value of Mh 2 (GKll) /Mh 2 (K13) seen at a fixed gas 
metallicity is large. Eow example, at Z fa 0.25 Zq, the median 
value of Mh2(GK11)/Mh 2(K13) is Ri 8, but the 16**' and 84*** 
percentiles are « 1.5 and ~ 38, respectively, revealing that there is 
another important parameter driving the differences seen in the H 2 
mass resulting from applying the GKl 1 and K13 prescriptions. 

To gain insight into what is driving the large dispersion 
in Mh 2 (GK11)/Mh 2 (K13) at a fixed gas metallicity, we se¬ 
lect galaxies in a narrow range of metallicities (shown by the 
hatched strip in the left-panel of Eig. lA2t and plot for those 
galaxies the ratio Mh 2 (GK11)/Mh 2 (K13) as a function of SER 
in the right-hand panel of Eig. IA2I It can be seen that galaxies 
with SFR< 0.01 Mq yr“^ exhibit the largest differences between 
MhjCGKH) and Mh 2(K13). We conclude that large deviations 
between H 2 mass calculated using the GKll or the K13 prescrip¬ 
tions appear in regimes of low metallicity and weak interstellar ra¬ 
diation fields. Note there is still a large scatter in the distribution 


Figure A2. The ratio between the H 2 mass calculated using the 
GKll prescription and the mass calculated using the K13 prescription, 
Mhj (GK11)/Mh2 (K13), as a function of the neutral gas metallicity 
at z = 0 in the simulation Ref-L100N1504 (left-hand panel). There is 
a trend of increasing Mh 2 (GK 11 )/Mh 2 {K13) with decreasing metal¬ 
licity. Lines with error bars show the median and the 16**' and 84**' per¬ 
centiles. We select galaxies in a narrow range in metallicity, shown by the 
hatched strip in the left-hand panel, and show in the right-hand panel the 
ratio ATh 2 (GK11 )/ATh2 (K13) as a function of the SFR. For reference, 
we show in both panels the ratio Mjjj (GK 11 )/Mh 2 (K13) = 1 as hor¬ 
izontal line. The major driver of the dispersion at a fixed metallicity is the 
SFR. 


shown at the right panel of Fig. lA2l which is due to the SER enter¬ 
ing in the H 2 fraction equations as SFR surface density. 

At 2 > 1 an interesting regime appears, which is the one of in¬ 
tense interstellar radiation fields, Gq > 100 in units of the Habing 
field (see Eig. ll2t . In this regime we observe that the KI3 prescrip¬ 
tion gives similar or slightly higher H 2 fractions than the GKll. 
Fig. lA3l is analogous to Fig. lA2l but for 2 « 2. As can be seen from 
the right-hand panel of Fig. lA3l galaxies with SFR> 1 Mq yr“*, 
have a ratio Mh 2 (GK11)/Mh 2 (K13) ~ 1. This is the regime of 
intense interstellar radiation fields. Note that for the same metallic¬ 
ity, galaxies at 2 = 0 have a higher ratio Mh 2 (GK11)/Mh 2 (K13) 
than galaxies at 2 = 2, on average. 

A3.1 Comparison with previous work 

In previous wo r ks, the empirical prescription of 
iBlitz & Rosolows^ (I 2 OO 6 I) has be en applied to gas particles 
in hydrodynamic simulations (e.g. iDuffv et akl l2012bl) . but re- 
stricting th e calcul ation to gas particles that are star-forming. 
IPuffv et^ (l2012bh adopted the following pres cription (following 
the observations presented in lLerov et alj200^ . 

where kh is Boltzmann’s constant. To get a better insight into the 
main differences between the theoretical prescriptions adopted here 
and the empirical prescription adopted in Duffy et al., we first esti¬ 
mate how much of the H 2 in the Ref-L100N1504 simulation with 
the prescriptions GKll and K13 applied is locked up in particles 
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Figure A3. As Fig. IA2I but for galaxies at 2 : = 2. Note that to produce 
the right-hand panel we select galaxies with lower metallicities than in the 
right-hand panel of Fig. lA2l to increase the number of galaxies selected. 


that are star-forming. We find that on average 82% of the H 2 is 
located in gas particles that are star-forming when the GKll pre¬ 
scription is applied. This percentage rises to 99% when the K13 
prescription is applied. However, we find that variations around that 
median are large, and galaxies can have as little at 20% of their H 2 
mass in gas particles that are star-forming. 

Second, we study the dependency of /hj, calculated using 
the GKll and K13 prescriptions, on pressure in the simulation 
Ref-L100N1504 and compare it with Eq. IA18I This is shown in 
Fig. lA4l The relations resulting from applying the GKll and K13 
prescriptions deviate significantly from Eg. lAlSI narticularlv in re¬ 
gions of low metallicity. This is a good example of the differences 
seen between applying the empirical relation of Eg. lAlSI and more 
physically-motivated prescriptions that include a dependence on 
gas metallicity. 


APPENDIX B: STRONG AND WEAK CONVERGENCE 
TESTS 

Here we compare the simulations Ref-L025N0376, Ref- 
L025N0752 and Recal-L025N0752 to address the effect of 
resolution on the H 2 masses of galaxies in a fixed cosmological 
volume. The comparison at a fixed volume is needed as simulations 
with different cosmological volumes sample different halos, which 
has the problem that differences can come from either resolution 
or halo sampling. We remind the reader that the simulation 
Recal-L025N0752 has four parameters that are modified slightly 
compared to Ref-L025N0376 and Ref-L025N0752 (see § for 
details). 

Eig. lBlI shows the H 2 mass as a function of total subhalo mass 
(dark plus baryonic) at 2 = 0 for the simulations Ref-L025N0376 
(solid lines), Ref-L025N0752 (dotted lines) and Recal-L025N0752 
(dashed line) using either the GKll (top panel) or K13 (bottom 
panel) prescriptions. The H 2 mass is measured in an aperture of 
30 pkpc, while Mxot is the total subhalo mass. We decide to use 
Mxot to characterise when the resolution becomes important in the 
H 2 mass of galaxies in the intermediate-resolution simulations be- 
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Figure A4. The fraction as a function of pressure for the simulation 
Ref-L100N1504 when the GKll (top panel) and the K13 (bottom panel) 
prescriptions are applied. Pixels are coloured by the median gas metallicity, 
as labelled at the top. The solid line shows the relation of Eq. IaTsI 


cause the subhalo mass is a well converged quantity. Eor the strong 
convergence test we need to compare Ref-L025N0376 with Ref- 
L025N0752 (i.e. fixed subgrid parameters). We see that the simu¬ 
lation Ref-L025N0752 shows H 2 masses for a fixed Mxot that are 
a factor of 2 — 3 larger than the H 2 masses in Ref-L025N0376 for 
Afxot > 5 X 1 0^° M(7). This differen ce is similar to the difference 
seen by S15 and lEurlong et al.l ( l2014h in the z = 0.1 galactic stellar 
mass function in the same set of simulations. 

For the weak convergence test we compare Ref-L025N0376 
with Recal-L025N0752. We find that the medians are closer to 
each other than in the strong convergence test. The difference seen 
between the median relations in the simulations Ref-L025N0376 
and Recal-L025N0752 is a factor of < 1.5 — 2 for Mxot > 
5 X 10^° M 0 , which again is consistent with the differences seen i n 
the galactic stellar mass function in S15 and lEurlong et al.l (|2014|) . 
Note that the exact prescription we use to calculate H 2 makes lit¬ 
tle difference to the offset obtained between the simulations we are 
comparing here. 

The higher resolution of Ref-L025N0752 and Recal- 
L025N0752 enables us to establish that, at the resolutions of the 
simulations Ref-L025N0376 and Ref-L100N1504, we can trust the 
H 2 content of subhalos with Mxot > 10^° M©. In the three sim¬ 
ulations analysed here there are ~ 1500 subhalos with Mxot > 
10^° Mq. 

Note that it is difficult to establish a H 2 mass below which res¬ 
olution effects become significant. This is because each gas particle 
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Figure Bl. The H 2 mass as a function of the total matter, Mxot (dark 
plus baryonic matter), at 2 = 0 for the simulations Ref-L025N0376, Ref- 
L025N0752 and Recal-L025N0752, as labelled, using the GKl 1 (top panel) 
and K13 (bottom panel) prescriptions. The H 2 mass here corresponds to the 
mass measured in an aperture of SOpkpc, while Mxot is the total subhalo 
mass. Lines show the medians, while the shaded region and errorbai's show 
the uncertainty on the median (e.g. the 16^^ and percentiles of the dis¬ 
tributions divided by where N is the number of objects at each bin). 

The simulations Ref-L025N0376 and Recal-L025N0752 agree very well at 
Mxot > 10^^ Mq regardless of the prescription used to calculate H 2 . At 
^Tot < 10^*^ M 0 the simulation Ref-L025N0376 starts to differ signifi¬ 
cantly from the higher resolution simulations which implies that resolution 
problems become important. 


has a weight, the H 2 fraction, and therefore a low H 2 mass can be 
obtained by either having a small H 2 fraction but a large number of 
gas particles, or by having a larger H 2 fraction but few gas particles. 
Thus, we find it more appropriate to express the resolution limit in 
terms of Mxot • We adopt Mxot = 10^*^ Mq as our resolution limit, 
and show unresolved results in Fig. [3] as a dotted line. In terms of 
H 2 mass, this resolution limit corresponds to ^ 5 x 10^ M© 
on average for the intermediate-resolution simulations, both for the 
GKll and K13 prescriptions. 
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